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SECTION  I 
INTRODUCTION 


2n  recent  years  a  great  deal  of  effort  has  gone  into  the 
detailed  computation  of  the  structure  of  atmospheric  nuclear  bursts. 

The  numerical  procedures  developed  for  these  calculations  assume  the 
existence  of  conventional  hydrodynamic  collieional  coupling  of  the  explo¬ 
sion  debris  to  the  ambient  atmosphere,  thus,  strong  shock  waves  are 
assumed.  This  assumption  certainly  is  valid  for  altitudes  below  about 

80  kilometers  or  so.  and  the  numerical  results  for  low  altitude  explosions 
appear  valid. 

The  assumption  of  strong  collieional  coupling  of  the 
debris  to  the  ambient  atmosphere  is  not  necessarily  valid  at  higher 
altitudes,  however,  as  a  result,  present  computational  techniques  do 
not  give  reliable  information  concerning  the  effects  of  high  yield,  high 
altitude  explosions. 


The  information  which  is  required  and  which  existing 
techniques  do  not  give  concerns  high  yield  explosions  at  high  altitudes. 
The  means  by  which  such  explosions  couple  to  the  atmosphere  and  the 
gross  effects  of  the  explosions  are  to  be  determined.  In  particular, 
techniques  capable  of  computing  the  transient  structure  of  the  expanding 
debris  of  an  explosion  in  the  megaton  range  at  an  altitude  between  100 
and  500  km  are  desired.  The  computed  structure  is  to  include  late-time 
effects,  including  disturbances  of  the  earth’s  magnetic  field. 


2.  REVIEW  OF  EXISTING  METHODS 


Existing  attempts  to  treat  high-altitude  explosions  either 
assume  conventional  coupling  and  apply  ordinary  hydrodynamic  computer 
codes  with  a  more  rarefied  atmosphere,  or  assume  there  is  no 
coupling  until  the  debris  has  expanded  to  some  "coupling  length,  "  after 
which  conventional  hydrodynamic  coupling  is  invoked.  The  assumption 
that  conventional  collisional  coupling  always  exists  obviously  breaks 
down  at  some  point:  at  a  sufficiently  high  altitude,  the  mean  free  path 
for  ambient  ions  becomes  greater  than  the  debris  radius,  and  the  con¬ 
ventional  momentum  transfer  process  disappears.  The  range  of  altitudes 
in  which  this  coupling  breakdown  occurs  is  not  known;  however,  because 
the  energy  released  by  the  explosion  so  disturbs  the  ambient  temperature 
and  ionization  that  estimates  of  the  effective  moan  free  path  for  momentum 
transfer  are  not  reliable. 

Existing  methods  of  accounting  for  the  coupling  break¬ 
down  do  so  by  attempting  to  estimate  an  effective  mean  free  path  or 
'coupling  length.  "  Various  processes  have  been  postulated  as  ultimate 
sources  of  the  coupling:  examples  include  the  "Longmire  Piston,  "  in 
which  the  motion  of  the  ionized  debris  generates  a  magnetic  field  that 
turns  or  "picks  up  "  the  ambient  ions;  two- stream  instability,  in  which 
the  electric  field  between  the  debris  ions  and  the  ambient  ions  grows  in 
an  unstable  fashion  until  it  becomes  strong  enough  to  pick  up  the  ambient 
ions,  and  the  collisional  model,  in  which  the  mean  free  path  for  direct 
collisions  is  estimated  on  the  basis  of  the  relative  velocities  of  the  debris 
and  ambient  ions.  All  of  these  postulated  coupling  processes  require 
rather  involved  mathematical  formulations  which  are  not  usable  in 
existing  hydrodynamics  codes.  As  a  result,  the  only  calculations  which 
have  been  made  with  the.se  processes  are  approximate  analytical  calcula¬ 
tions  solely  intended  to  demonstrate  the  possibility  that  the  process  is 


significant.  No  detailed  calculations  intended  to  show  the  actual  effect 
of  such  processes  have  been  made. 

3*  SUMMARY  of  the  present  approach 

Because  of  the  potential  significance  of  coupling  pro¬ 
cesses  of  the  type  discussed  above,  a  physics  code  is  developed  which 
can  handle  the  mathematical  formalism  associated  with  these  coupling 
processes.  A  suitable  coding  technique  exists  and  has  been  used  on 
transient  interaction  problems  in  the  past;  further,  the'code  has  been 
redeveloped  to  the  point  where  it  is  useful  for  studying  the  detailed 
behavior  of  the  proposed  coupling  models.  This  code,  which  is  an 
Eulerian  code  with  a  floating  mesh,  employs  the  characteristic  flux 
differencing  technique:  the  changes  in  all  properties  during  a  time- 
step  are  caused  by  the  fluxes  of  those  properties  across  the  surface  of 
the  cell  during  the  time-step.  Because  surface  fluxes  in  a  rarefied  gas 
do  not  involve  thermodynamics  or  kinetic  processes  such  as  collisions, 
it  is  possible  to  construct  the  entire  finite -difference  code  without  such 
models.  These  models  are  added  in  subroutine  form  for  the  evaluation 
of  properties  such  as  partial  pressures  and  sound  speeds  within  each 
cell.  This  differencing  procedure  has  the  great  versatility  of  permit¬ 
ting  a  complete  change  of  the  kinetic  and  thermodynamic  models  to  be 
made  merely  by  changing  a  few  subroutines. 

This  code  structure  will  be  useful  for  carrying  out 
detailed  one-dimensional  studies  of  the  interaction  of  debris  and 
ambient  ions  for  each  of  the  postulated  coupling  processes.  The  results 
of  the  more  effective  of  these  processes  can  be  scaled  into  "coupling 
laws,"  which  might  be  used  in  a  two-dimensional  calculation  of  the  gross 
expansion.  In  addition,  two-dimensional  studies  can  be  made  for  those 
coupling  mechanisms  which  are  inherently  two-dimensional;  the  "Longmire 
Piston"  with  its  streaming  debiis,  perpendicular  magnetic  field,  and 
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turning  of  the  ambient  ions  is  an  example  of  a  coupling  process  which 
must  be  treated  as  at  least  two-dimensional. 

In  addition,  numerical  studies  for  evaluating  the  several 
coupling  mechanisms  and  calculating  the  gross  expansion  structure  of 
high-altitude  explosions  can  l?e  carried  out. 


SECTION  II 


FORMULATION  OF  THE  PROBLEM 
AND  METHOD  OF  SOLUTION 


The  problem  being  treated  here  is  that  of  developing  a 
computer  code  which  can  be  used  for  the  dual  purposes  of  (1)  testing 
the  efficacy  of  proposed  mechanisms  for  coupling  the  debris  from  a 
high-altitude  explosion  to  the  ambient  atmosphere,  and  (2)  computing 
the  overall  features  of  explosions  for  which  the  coupling  mechanisms 
have  been  determined.  The  method  of  treating  this  problem  consists 
of  using  the  governing  equations  in  a  form  that  is  valid  regardless 
of  the  coupling  mechanism,  selecting  an  appropriate  numerical  pro¬ 
cedure  for  solving  these  equations,  and  developing  (and  coding  as  sub¬ 
routines)  models  of  the  proposed  coupling  mechanisms. 


A  mathematical  description  of  explosion  phenomenology 
must  involve  conservation  of  mass,  momentum,  and  energy;  probably 
must  include  radiative  phenomena;  and  must  account  for  electromagnetic 
effects.  In  addition,  a  set  of  constitutive  equations  is  required,  includ¬ 
ing  equations  of  state  and  of  electrical  conductivity.  Finally,  appropriate 
boundary  conditions  must  be  prescribed. 


The  interacting  explosion  debris  and  ambient  atmosphere 
may  consist  of  molecular,  atomic,  and  ionized  species  plus  electrons. 

It  is  assumed  that  there  are  N  such  species  present,  including  elec¬ 
trons;  an  individual  species  is  denoted  by  subscript  alpha.  The  symbols 


used  in  the  following  equations  are  defined  in  the  Glossary,  and  the 
normalization  is  discussed  in  Subsection  II-2  below. 


a*  Hydrodynamic  Conservation  Equations 

The  equations  expressing  the  conservation  of  mass, 
momentum,  and  energy  of  species  a  are  obtained  by  taking  moments 
of  the  Boltzmann  equation  for  species  a  (Appendix  I).  The  resulting 
conservation  equations  are 


Mass: 


ap, 


a  - 
+  v 


at 


Momentum: 


Kma»a(1> 


^(Pa<“a>)+''(Pa<“a><“a>+Pa) 

=  Pa<  ?a /ma>  +  K  wa(u) 

Energy: 

£(Paea)  +  *  •  (Pa£a<%>  +  Pa  *  <^>) 

=  -  ^  *  Qa  +  Pa( Fa  /ma)-  <^> 

+  Kmawa(e  +  u2/2) 


(1) 


(2) 


(3) 


Equations  1  through  3  are  valid  for  nonequilibrium 
flow  as  well  as  equilibrium  flow.  The  wa  terms  on  the  right-hand 
sides  of  these  equations  represent  sources  and  sinks  of  mass,  momen¬ 
tum,  and  energy  resulting  from  particle  collisions.  In  the  event  of 
equilibrium  flow,  these  source  terms,  which  derive  from  the  collisional 
term  in  the  Boltzmann  equation,  vanish,  and  the  pressure  tensor 
is  given  by  the  equilibrium  equation  of  state.  In  the  event  of 
nonequilibrium  or  of  chemically- reacting  flow,  these  terms  are 
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evaluated  by  approximating  the  colliaional  terms  in  the  Boltzmann 

equation  by  experimentally  determined  eross  sections  for  the  collision 
in  question  (see  Appendix  I). 


Radiation  Transport  Equations 


Equations  for  the  energy  density  and  heat-flux  vector 
of  radiation  of  frequency  v  are  obtained  by  taking  moments  of  the 
Boltzmann  equation  for  the  distribution  function  of  photons  of  energy 
hv  (Appendix  III).  The  resulting  equations  are 


1  9eRV  -a  - 

- rr~  +  v  •  qr 

c  at  Rv 


1 

c 


at 


+  v  • 


-  Bu?*>a[’%eRv(1  ’ Jva)-jvo]  (4) 

-  Bu  Qj^  §Pa HVa (l  -  JVa  +  SVa]  (5) 


Equations  4  and  5  are  valid  for  any  radiation  con¬ 
dition.  The  terms  on  the  right-hand  sides  of  these  equations  represent 
sources  and  sinks  of  radiative  energy  and  heat  flux  resulting  from 
absorption,  emission,  and  scattering. 


Maxwell's  Equations 


Maxwell's  equations  must  be  normalized  into  a  form 
which  can  be  related  to  a  wide  variety  of  initial  conditions  in  order  to 
be  convenient  for  general  computational  purposes.  It  is  shown  in 
Appendix  IV  that  either  the  rationalized  mks  form  or  the  Gaussian 
form  of  Maxwell's  equations  can  be  normalized  into  the  following  set: 
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V  x  E 


(8) 


V  X  B 


1  SB 

c  at 

£j  l  9E 

c  cat 


(9) 


Equations  6  through  9  are  valid  for  the  electro¬ 
magnetic  fields  in  any  moving  medium  with  no  inaccessible  charges 
or  currents;  that  is,  all  charges  and  currents  are  assumed  accounted 

for  by  the  moving  particles  of  the  medium,  and  there  are  no  polariza¬ 
tion  currents. 


Constitutive  Equations 

The  field  equations  listed  above  must  be  complemented 
by  a  set  of  constitutive  equations  which  govern  the  transport  properties 
of  the  fields  in  the  medium  of  propagation.  Thus,  Equation  2 
requires  an  equation  for  the  partial  pressure  tensor  of  species  a.  (or 
for  the  transport  of  momentum)  and  Equation  3  requires  an  equation 
for  the  heat-flux  vector  of  species  a  (or  for  the  transport  of  energy). 
Similarly,  Equation  5  requires  an  equation  for  the  radiative  pres¬ 
sure  tensor  (or  for  the  transport  of  radiative  momentum),  while 
Equation  9  requires  an  equation  for  the  current  (or  for  the  transport 
of  charge).  These  constitutives  represent  moments  of  the  nonequi¬ 
librium  distribution  functions,  and  generally  are  not  known.  However, 
it  is  possible  to  postulate  various  models  for  these  needed  equations, 
and  it  will  be  shown  (Subsection  IH-2c)  that  these  constitutive  equations 
do  not  affect  or  enter  into  the  differencing  technique  for  the  numerical 
code.  As  a  result,  it  is  possible  to  change  the  constitutive  equations 
used  in  the  code  merely  by  changing  the  corresponding  subroutines. 

A  number  of  models  for  the  hydrodynamic  constitutive 
equations  are  available;  these  models  usually  are  based  either  on  the 
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apman  and  Enskog  perturbation  theory  for  the  nonequilibrium  die- 
tribution  function  (Reference  1).  or  on  extensive  experimental  data. 
The  simples,  models  having  any  physical  reality  for  the  hydrodynamic 

— r  'T,i0nS  <Md  the  —  '«  *»•  initial 

numexical  work  with  the  code)  are 

P“*  =  Xk*Ta? 


Equation  lo  is  valid  equation  for  nonequilibrium  flow,  but  is  an 
approximate  model  because  shearing  stresses  have  been  ignored 

(see  for  example.  Section  2.  Chapter  IX  of  Reference  2).  normal- 
ized  form.  Equations  10  and  11  are 
=» 

-  P  =  n  KT  I 

a  a  a  /i?i 


......  Q„ 

a  . 


A  satisfactory  constitutive  equation  for  the  radiative 
pressure  tensor  appears  to  be  that  given  by  the  Milne-Eddington 
approximation  (which  has  been  used  successfully  in  astrophysics,-. 

4  =  2*? 

3  1  (1 

«  can  be  shown  that  Equations  4  and  5  with  Equation  14  lead 
to  the  correct  radiation  transport  formulae  in  the  optically  thick  and 
optically  thin  limits  as  well  as  in  the  limit  of  isotropic  radiation. 

quations  4  and  5  with  Equation  14  are  similar  in  form  to  the 
equations  obtained  in  the  firs,  approximation  in  the  spherical- 
harmonic  method  of  neutron  transport  theory;  in  that  theory  it  is 
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known  that  the  odd  approximations  are  more  accurate  than  the  suc¬ 
ceeding  even  approximations  (Reference  2).  A  final  comment 
on  the  validity  of  Equation  14  is  that  it  deals  with  the  transport  of 
radiative  momentum,  and  that  radiative  momentum  considerations 
only  enter  in  a  relativistic  theory.  Consequently,  Equation  14 

is  expected  to  be  a  satisfactory  constitutive  equation  in  a  nonrelativistic 
theory. 


The  necessary  constitutive  equation  for  the  current 
density  in  Equation  9  is  straightforward  as  long  as  all  of  the  charges 
move  by  convective  flow  and  no  conduction  currents  exist;  this  equation 
is 

>  =  £W3a>  (15| 

If  problems  involving  conduction  currents  are  to  be  treated,  it  will  be 
necessary  to  develop  an  expression  for  the  conductivity  oi  the  flowing 
gas.  The  initial  numerical  work,  however,  will  rely  on  Equation  15. 


e*  Coupled  Equations 

The  momentum  and  heat  transfer  effects  governed  by 
the  radiation  transport  equations  and  Maxwell’s  equations  must  be 
coupled  into  the  hydrodynamic  equations  to  provide  a  complete  des¬ 
cription  of  the  physics  of  the  problem. 

As  shown  in  Appendix  IV,  the  electromagnetic  com¬ 
ponent  of  the  force  term  in  Equation  2  has  the  normalized  form 


(16)' 


normalized  work  term 


energy  is  given 


Similarly,  the  electromagnetic  component  of  the 
in  Equation  3  is  shown  to  be 

v  -  -  _ 

Vua  = 


The  heat  added  to  species  a  by  absorption  of  radiative 
by  Equation  273  of  Appendix  III  as 

*  Bo Bu/[eRv’tRvf‘  -  JVa) - Jvjdv 


Coupling  Equations  16  through  18  into  Equations  1 
through  3  gives  the  coupled  hydrodynamic  equations: 

bf> 

It  +M',a<3a>)  *  K%»a(l)  (1 

S('’a<%>)  +  ?-(pa(^X5a>  +  4) 

=  o  ii  i  %  /-  <3a)xB\l 

a|ma  *!  <»a  VE+  c  /J+  Kma"a(u)  <2( 

^(paea)  +  ?  '  (0oea<3a>  +  Pa-  <ua>) 


■  v  •  +  bu/[=Rv  "R^l  -  Jvj  -  jv Jdv 

JJT  <ua>  (?a  +  -jjf- )  +Km0w0(e  +  u2/2) 


The  force  Fa  in  Equations  20  and  21  ls  included  to 

prov.de  an  easy  means  of  adding  the  gravitational  acceleration  a,  a 
future  date. 
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f. 


Boundary  Conditions 


Appropriate  boundary  conditions  must  be  provided  on 
the  boundary  of  the  region  of  interest.  In  mathematical  terms,  the 
boundary  requirement  is  that  all  characteristics  crossing  into  the 
region  of  interest  from  the  exterior  region  must  be  known.  In  the 
case  of  Maxwell's  equations,  this  statement  means  that  either  there 
must  be  no  charges  or  currents  in  the  exterior  region  (in  which  case 
the  boundary  conditions  are  those  of  the  fields  at  infinity)  or  the  time- 
dependent  fields  on  the  boundary  must  be  prescribed  ahead  of  time. 

The  boundary  conditions  for  the  radiation  transport  equations  are 
similar  to  those  for  Maxwell's  equations:  either  there  must  be  no 
emission,  absorption,  or  scattering  in  the  exterior  region  (in  which 
case  the  boundary  condition  is  that  of  constant  flux  from  infinity)  or 
the  time- dependent  radiation  fields  on  the  boundary  must  be  prescribed 
ahead  of  time. 


The  boundary  conditions  required  for  the  hydrodynamic 
equations  depend  upon  whether  the  flow  within  the  boundary  is  super¬ 
sonic  or  subsonic  relative  to  the  boundary.  If  the  flow  is  supersonic 
relative  to  the  boundary  and  crossing  out  of  the  region  of  interest,  there 
can  be  no  characteristics  entering  the  region  of  interest  (hydrodynamic 
characteristics  propagate  at  the  speed  of  sound)  and  no  boundary  con¬ 
ditions  are  needed.  In  any  other  case,  however,  the  flow  on  the 
exterior  of  the  boundary  must  be  prescribed. 

For  all  of  the  numerical  calculations  planned  for  the 
near  future,  these  boundary  conditions  will  be  satisfied  by  assuming 

that  the  fields  at  the  boundary  are  identical  to  the  undisturbed  fields 
at  infinity. 
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l'  NORMALIZATION 


The  normalization  of  the  equations  listed  in  Subsection 
1-1  is  indicated  in  the  notations:  tnis  normalization  is  intended 

to  ensure  that  the  only  quantities  in  Equations  4  through  9  and  12 

through  21  u;hich  are  not  of  order  unity  are  the  parameters  c,  B 

K.  K  B„,  and  K.  These  parameters  are  measures  of  the  rela¬ 
tive  importances  of  the  various  physical  processes  represented  in 
these  equations.  Thus,  c  is  the  ratio  of  the  speed  of  light  to  the 
characteristic  hydrodynamic  velocity,  and  is  a  measure  of  the  impor¬ 
tance  of, .transient  radiative  and  electromagnetic  effects.  Large  c 
implies  fhat  radiative  energy  transport  occurs  on  a  time  scale  which 

13  Sh°rt  C°mpared  the  hydrodynamic  time  scale,  and  therefore.these 
effects  become  quasi-steady,  .mall  c  (c  of  order  unity,  implies  that 
radiative  and  electromagnetic  transient  effects  are  important,  and 
therefore,  the  radiative  and  electromagnetic  energy  densities  in 
transit  are  important.  Bu.  the  Bouguer  number,  is  a  measure  of  , he 
relative  opacity  of  the  radiating  gas.  Large  Bu  implies  there 
are  many  optical  path  lengths  in  a  characteristic  geometric  length, 
and  therefore, absorptive  and  emissive  properties  of  the  gas  are 
important.  Small  Bu  implies  that  a  characteristic  geometric  is  short 
compared  to  an  optical  path  length,  and  therefore.the  absorptive  and 
emissive  properties  of  the  gas  are  unimportant,  g .  the  ratio  of  the 
geometric  scale  length  to  the  Debye  length,  is  a  measure  of  the  rela¬ 
tive  lengths  of  the  hydrodynamic  and  electromagnetic  fluctuations. 

Large  X  implies  that  a  geometric  scale  length  contains  many  fluctu¬ 
ations  of  the  electromagnetic  fields,  and  therefore.the  charge  and 
current  densities  (the  sources  of  the  fluctuations)  are  important. 

Small  1  implies  that  the  electromagnetic  fluctuations  are  not  impor- 
tant  over  a  geometric  scale  length,  and  therefore,  the  charge  and 


current  densities  are  unimportant.  X,  which  is  very  nearly  the 
square  of  th> !  reciprocal  of  the  Mach  number,  is  a  measure  of  the 
relative  energies  of  random  motion  and  directed  motion  in  the  gas. 
Large  K  implies  that  the  thermal  energy  of  random  motion  is  large 
compared  to  the  directed  energy  of  the  streaming  motion,  and  there¬ 
fore,  pressure  effects  are  important  Small  X  implies  that  the 
thermal  energy  of  the  random  motion  is  small  compared  to  the 
directed  energy  of  the  streaming  motion,  and  therefore,  j  -essure 
effects  are  unimportant.  7t[  is  a  measure  of  the  relative  importance 
of  the  hydrodynamic  momentum  of  ,a  particle  and  of  the  impulse 
delivered  to  the  particle  in  unit  time  by  the  local  electric  field.  Thus, 
^■ar6e  ^  implies  that  electromagnetic  effects  are  not  important, 
while  small  fl{  implies  that  such  effects  become  dominant.  B0,  the 
Boltzmann  number,  is  a  measure  of  the  relative  importance  of  radiative 
energy  transport  and  of  convective  energy  transport.  Large  BQ  signi¬ 
fies  that  radiative  transport  is  important,  while  small  BQ  signifies 
that  radiative  transport  is  unimportant.  K  is  a  measure  of  the  rela¬ 
tive  importance  of  the  change  in  particle  density  of  a  species  caused 
by  chemical  reactions  and  the  change  caused  by  hydrodynamic  convec¬ 
tion.  Chemical  reactions  are  important  for  large  K  and  unimportant 
for  small  K. 

3. 

The  problem  outlined  in  Subsection  II-l  is  of  such  com¬ 
plexity  that  only  a  numerical  solution  appears  possible.  There  are  a 
number  of  finite- difference  methods  available,  but  nearly  all  of  these 
methods  have  limitations  which  preclude  the  possibility  of  treating 
the  above  equations  in  their  full  generality.  Such  limiting  techniques 
include  the  conventional  Eulerian  and  Lagrangian  procedures  as  well 
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as  the  Particle-in- Cell  approach.  By  reviewing  these  procedures, 
however,  it  is  possible  to  define  the  desirable  features  of  a  more  ” 
general  numerical  approach,  and  then  to  outline  such  an  approach. 

a*  Eulerian  MethnH« 

Eulerian  finite -difference  methods,  which  are  frequently 
used  for  two-dimensional  problems,  proceed  by  applying  the  governing 
1  ferential  equations  at  a  set  of  discrete  points  in  space.  The  differ- 
entials  in  the  equations  are  approximated  by  finite-differences  between 
the  fixed  points.  Because  of  the  fixed  geometry  associated  with  these 
methods.  Eu'.erian  procedures  are  readily  able  to  compute  shearing 
stresses  ani  to  treat  tangential  forces.  This  capability  is  a  necessity 
for  the  complex  problems  considered  here  because  of  the  shearing 
stresses  associated  with  magnetic  field  and  the  resulting  nonequality 
of  the  diagonal  components  of  the  hydrodynamic  pressure  tensor. 

One  of  the  major  disadvantages  associated  with  Eule.  ian 
ethods  is  that  they  are  unable  to  maintain  sharp  contact  discontinuities 
between  two  separate  fluids.  This  inability  arises  because  the  fluid 
characteristics  at  a  point  represent  the  average  characteristics  of  the 
fluid  in  a  small  volume  surrounding  the  point;  thds.  when  a  contact 
discontinuity  enters  such  a  volume,  the  sharp  discontinuity  is  repre¬ 
sented  by  an  average  a,  the  point  a,  the  center  t^f  the  volume,  m 
this  fashion,  the  contact  discontinuity  becomes  smeared  out  after  a 
few  computational  cycles.  A  detailed  study,  including  a  numerical 
example,  of  this  smearing  of  contact  surfaces  is  given  in  Subsection 
Of  Appendix  F  of  Reference  3.  Unfortunately,  many  of  the  problems 
considered  here  have  interacting  gas  streams  or  colliding  shock  waves 
w  ich  produce  contact  surfaces;  consequently,  this  inability  of  Eulerian 

methods  to  properly  represent  such  contact  surfaces  is  a  strong  disad- 
vantage  of  these  methods. 


A  second  disadvantage  associated  with  Eulerian  methods 
is  that  they  generally  are  unable  to  treat  non- rectangular  boundaries. 
This  inability  results  from  the  finite-difference  approach  of  representing 
differentials  by  differences  taken  between  successive  points  of  the  com¬ 
putational  net.  It  is  straightforward  to  treat  tl  e  problem  of  rectangular 
walls  by  imposing  symmetry  on  the  fixed  points  in  the  computational 
network,  but  similarly  satisfying  procedures  for  other  geometries  are 
lacking.  In  the  problems  considered  here,  however,  the  constraints 
impos 3d  by  the  atmospheric  variations  and  the  earth's  magnetic  field 
may  dictate  a  nonrectangular  boundary  for  the  numerical  pr  oblem. 

As  a  result,  this  limitation  of  Eulerian  methods  is  a  disadvantage. 

b*  Lagrangian  Methods 

Lagrangian  computational  methods  refer  to  a  specific 
set  of  fluid  particles  rather  than  to  a  set  of  spatial  points.  Newton's 
laws  are  applied  to  these  representative  points  to  obtain  the 
motion  of  the  fluid  as  a  whole.  Because  i..?grangian  methods  follow 
individual  fluid  particles,  such  techniques  are  able  to  represent  con¬ 
tact  surfaces  between  two  gases;  as  discussed  above,  this  ability  is 
important  for  the  problems  to  be  treated. 

The  major  disadvantage  of  Lagrangian  methods  is  that 
they  do  not  treat  shearing  stresses  satisfactorily.  This  difficulty 
arises  because  the  calculation  of  the  stresses  is  based  upon  the  geo¬ 
metric  distances  between  the  fluid  particles,  and  the  calculation  of 
these  distances  becomes  unacceptably  tedious  after  the  fluid  has  under¬ 
gone  some  distortion  because  of  the  shearing  stresses.  As  discussed  in 
Subsection  II- 3a,  however,  shearing  stresses  are  important  in  the 
present  application,  and  the  inability  of  Lagrangian  methods  to  handle 
them  is  a  serious  disadvantage. 
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The  Particle-in-Cell  (or  PIC)  technique  represents  an 

attempt  to  gain  the  advantage  both  of  Eulerian  and  Lagrangian  method, 

by  superimposing  a  set  of  Lagrangian  mass  points  on  a  set  of  Euieriar 

volume  elements.  The  fluid  stresses  are  computed  on  the  basis  of  the 

average  properties  of  all  of  the  mass  point,  in  each  of  the  volume 

elements,  and  the  motion  of  the  fluid  is  computed  by  using  these 

stresses  to  determine  the  acceleration  and  velocity  of  each  of  the 

mass  points^.  The  use  of  volume  elements  permits  this  technique  to 

treat  shearing  stresses,  while  the  presence  of  mass  points  allows 

good  representation  of  contact  surfaces.  As  a  result,  the  PIC  tech- 

mque  has  advantages  over  the  Eulerian  and  Lagrangian  methods  des- 
cribed  above. 


The  disadvantage  of  being  unable  to  treat  nonrectangule 
boundaries  remains,  however.  The  presence  of  this  disadvantage  is 
because  of  the  use  of  Eulerian  difference  procedures  for  computing  th, 
stresses  on  the  mass  points,  and  results  in  exactly  the  same  way  as 
the  rectangular  boundary  requirement  discussed  in  Subsection  U-3a. 

A  further  disadvantage  of  the  PIC  technique  is  that  the 
mass  points  are  par,  of  a  Lagrangian  net  and  there  may  be  only  a  few 
pomts  of  each  species  in  a  given  Eulerian  cell.  As  a  result,  it  is 
difficult  to  represent  the  smooth  change  in  the  density  of  a  species  bee, 
of  chemical  reactions  by  changing  the  number  of  mass  points  in  the 
cell.  Consequently,  the  PIC  technique  shares  the  Lagrangian  disadvan, 
age  of  not  being  well  suited  for  chemical  nonequilibrium. 

d*  .Characteristic  Flny  IvTethnH 

A  method  which  uses  moving  Eulerian  cell  boundaries 

to  follow  contact  discontinuities  is  aMn  a _ 

continuities  is  able  to  overcome  the  disadvantages 


of  all  Of  the  above  methods.  Such  a  method  is  discussed  in  great 
detail  in  Reference  3.  The  basic  concept  of  this  method  is  that  of 
choosing  an  Eulerian  spatial  mesh  and  applying  the  governing  equa¬ 
tions  in  integral  form,  rather  than  differential  form,  to  the  mesh. 

The  resulting  computational  equations  give  the  changes  in  fluid  proper¬ 
ties  within  a  cell  in  terms  of  the  fluxes  of  mass,  momentum,  and 
energy  of  the  fluid  crossing  each  face  of  the  cell.  Because  the  fluid 
properties  are  assumed  constant  within  each  cell,  however,  the 
properties  on  any  face  of  the  cell  must  be  given  in  terms  of  a  simple 
wave  emanating  from  within  the  cell  (Section  29,  Reference  4).  Thus, 
the  fluxes  crossing  the  cell  surface  are  computed  by  the  method  of 
characteristics,  and  no  distances  to  cell  centroids  are  needed. 
Furthermore,  the  determination  by  the  method  of  characteristics  of 
the  fluxes  crossing  each  cell  face  is  based  upon  velocities  relative  to 
the  cell  face;  thus,  the  Eulerian  grid  can  move  through  the  gas. 

A  computational  procedure  based  upon  this  character¬ 
istic  flux  method  has  all  of  the  advantages  of  the  above  methods  without 
any  of  the  disadvantages.  Thus,  the  use  of  integral  equations  replaces 
calculations  of  gradients  vith  calculations  of  forces  and  fluxes  on  cell  H 
surfaces.  Because  the  Eulerian  cell  surfaces  do  not  become  distorted  by 
the  shear  flow,  these  force  and  flux  calculations  remain  straightforward. 
At  the  same  time,  the  movable  cell  boundaries  permit  good  representa¬ 
tion  of  contact  discontinuities.  For  example,  cell  surfaces  can  be 
aligned  with  and  allowed  to  move  with  material  interfaces;  thereby, 
preventing  the  usual  smearing  of  the  interfaces  into  adjacent  cells. 

Such  a  technique  was  used  successfully  in  Reference  5  in  treating  the 
contact  surface  formed  by  the  colliding  shock  waves  when  a  reentry 
vehicle  penetrates  a  blast  wave. 


The  characteristic  flux  technique  also  is  able  to  treat 
chemical  nonequilibrium  without  difficulty.  Since  the  procedure  is 
an  Eulerian  procedure,  masses  are  computed  by  the  use  of  continuous 
densities  rather  than  by  the  use  of  discrete  mass  points.  As  a  result, 
reaction  rate  equations  can  be  coupled  into  the  procedure  directly. 

Finally,  the  characteristic  flux  procedure  can  treat 
boundaries  of  any  geometry.  This  ability  results  from  the  replace- 
ment  of  gradients  in  the  differential  equations  by  fluxes  crossing  cell 
walls  in  the  integral  equations.  Thus,  i,  only  is  necessary  to  align 
cell  surfaces  with  the  boundary  to  treat  arbitrary  geometries.  For 

example,  ellipsoidal  as  well  as  spherical  bodies  were  treated  in 
References  3  and  5. 


4*  METHOD  OF  SOLUrtriM 

As  discussed  in  Subsection  II- 3.  the  characteristic  flux 
numerical  method  is  the  only  available  procedure  which  can  treat  with 
the  required  generality  Equations  4  through  9  and  19  through  21 
along  with  the  associated  boundary  conditions.  Consequently,  these 
equations  are  integrated  over  a  cell  volume,  and  the  method  of  char¬ 
acteristics  is  used  to  obtain  the  required  fluxes  crossing  the  cell 
surfaces. 


The  Eulerian  cells  are  referred  to  a  moving  mesh. 

As  discussed  in  Subsection  11-3.  d,  such  a  mesh  makes  it  possible  to 
maintain  internal  interfaces  in  the  gas;  subroutines  are  provided  for 
determining  the  mesh  motion  required  for  the  following  of  the  inter¬ 
faces.  In  addition,  the  use  of  a  moving  mesh  makes  it  possible  to 
have  the  computational  mesh  expand  with  the  flow  (as  in  the  case  of 
the  expanding  debris  from  an  explosion);  thereby, permitting  a  constant 
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number  of  cells,  with  the  resulting  uniform  accuracy,  to  rover  the 
entire  flow  field.  This  type  of  motion  obviates  the  need  for  rezoning 
the  mesh  periodically  by  adding  or  deleting  cells. 

Finally,  all  of  the  calculations  involving  the  constitu¬ 
tive  equations  are  carried  out  in  subroutines;  thereby, facilitating  the 
changing  of  such  equations  as  new  models  are  developed. 


SECTION  III 


FORMULATION  OF  THE  PHYSICS  CODE 


A  general  physics  code  is  developed  to  treat  the  equa¬ 
tions  presented  in  Section  II  under  a  wide  selectipn  of  boundary  conditions. 
This  code  is  developed  by  integrating  the  governing  equations  over  the 
volume  of  a  cell  and  then  using  the  method  of  characteristics  for  the 
computation  of  the  field  properties  on  the  surfaces  of  the  cells.  The 
constitutive  equations  are  added  as  subroutines.  The  computational 
mesh  is  set  up  in  a  general  form  which  permits  selection  of  plane, 
cylindrical,  or  spherical  geometries  in  either  one-  or  two-dimensions 
as  needed.  A  method  for  causing  the  mesh  to  move  in  such  a  way  as 
to  follow  a  set  of  prescribable  interfaces  is  developed.  A  stability 

analysis  of  the  resulting  code  is  carried  out,  and  detailed  flow  charts 
of  the  code  are  presented. 


INTEGRATION  OF  THE  EQUATIONS 

The  equations  are  integrated  in  such  a  way  as  to  be 
correctly  represented  in  either  one-  or  two-dimensions  and  in  either 
plane,  cylindrical,  or  spherical  geometry.  These  three  coordinate 
systems  plus  a  typical  computational  cell  in  each  system  are  illustrated 
in  Figure  1.  The  geometrical  definitions  which  specify  a  cell  in  any  of 
these  systems  are  shown  in  Figure  2.  The  field  properties  represented 
by  Equations  4  through  9  and  19  through  21  are  assumed  to  be  con¬ 
stant  within  each  cell  during  a  time  step  and  constant  on  each  face  of 

each  cell  during  the  time  step.  The  equations  are  then  integrated  over 
the  volume  of  a  cell  and  over  a  time  step. 
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a*  Hydrodynamic  Conservation  Equations 

Equations  19  through  21  with  the  cell  integration 


indicated  are 


/dVX^+MVV)--*a(1>j  ■  0 
jH  dtp;Ma>) + ’ ♦  i)  -  ja(2j 


and 


XdVXd‘fe(Pae«)+  ;  '(wv  +  •  <V  +  5a) 
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where 
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V1) 

=  Kmawa(1) 
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3a<2> 

Kma*a(tt) 
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■V3> 

Bo  B„/ (‘ -  Jva)  -  Jvt 
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,  pa  ^  “4\ 

m  <V  \Fa  +TE/+Krna 

wa(e  +  u2/2) 

(27) 

The  integrations  in  Equations  22  through  24  are  carried  out  over  the 
volume  of  a  cell,  V,  and  over  a  time  step,  At.  In  carrying  out  these 
integrations,  it  is  assumed  that  the  time  rates  of  change  and  the  source 
functions,  Jr,  within  each  cell  are  constant  during  the  time  step.  The 


divergence  term  is  converted  to  an  integral  over  the  cell  surface  by 
means  of  Green's  theorem,  and  the  resulting  properties  on  the  cell 
surfaces  are  assumed  constant  during  the  time  step.  In  integrating 
the  unsteady  term  in  these  equations,  it  must  be  remembered  that 
the  cell  volume  is  a  function  of  time  because  of  the  moving  mesh;  thus, 
this  term  is  integrated  by  parts,  and  the  resulting  expression  for  the 
rate  of  change  of  cell  volume  is  replaced  by  an  integral  over  the  sur¬ 
face  of  the  cell  of  the  normal  component  of  the  velocity  of  the  surface. 
These  integration  procedures  are  carried  out  in  Appendix  V  for  repre¬ 
sentative  scalar  and  vector  equations. 


Making  use  of  the  model  integrals  provided  by  Equa¬ 
tions  3i9  and  331  and  by  Equations  320  ,  339  ,  and  369  of 

Appendix  V,  the  integrated  forms  of  Equations  22  through  24  become 
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and 
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In  Equations  28  through  31  ,  Wn  is  the  velocity  of  the 
cell  surface  along  its  outward  normal,  <ua>n  is  the  velocity  along  the 
outward  normal  of  the  cell  surface  of  species  a,  and  vr  and  vz  are 
the  direction  cosines  between  the  r  and  z  axes  and  the  outward  normal 
to  the  cell  surface.  The  index  (i+1/2,  j+1/2)  denotes  properties  evalu¬ 
ated  within  cell  (i+1/2,  j+1/2)  ;  raised  indices  denote  properties  evaluated 
at  the  end  of  the  time  step,  while  lowered  indices  denote  properties  evalu¬ 
ated  at  th^<)beginning  of  the  time  step. 

b*  Radiation  Transport  Equations 

Equations  4  and  5  with  the  cell  integration  indicated 

are 
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Equations  32  and  33  are  integrated  with  the  same 
assumptions  and  in  the  same  manner  as  Equations  22  through  24  .  The 
resulting  integrated  forms  are 
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c.  Maxwell's  Equations 


Equations  6  and  7  will  always  be  satisfied  if  they  ar 
satisfied  by  the  initial  conditions;  consequently,  they  need  not  be  differ 
enced.  Equations  8  and  9  with  the  cell  integration  indicated  are 


and 


/dv/tdt[lf +  c5x£]  =  0 

XdV-£dtM'c'x5+xT] -  ° 


Equations  39  and  40  are  integrated  with  the  same 
assumptions  and  in  the  same  manner  as  Equations  22  through  24  . 
The  resulting  integrated  forms  are 
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2‘  characteristic  solutions  ON  CELL  ^OUNDAT?!^ 

Equations  28  through  31  ,  36  through  38  ,  and  41 

through  44  all  contain  terms  evaluated  on  the  faces  of  the  computational 
cells.  These  terms  represent  fluxes  of  the  property  in  question  across 
the  cell  faces,  and  are  calculated  by  the  method  of  characteristics. 


a*  General  Theory 


The  various  field  properties  governed  by  the  above  inte¬ 
grated  equations  are  considered  to  be  constant  within  each  cell;  thus, 
the  cell  faces  constitute  a  region  bounding  a  region  of  constant  field.  As 
shown  in  Section  29  of  Reference  4,  however,  the  field  in  a  region 
adjacent  to  a  region  of  constant  state  must  be  a  simple  wave  emanating 
from  the  region  of  constant  state.  The  leading  edge  of  this  simple  wave 
propagates  along  the  undisturbed  characteristic  of  the  constant  state,  and 
in  the  limit  of  small  differences  between  th*  constant  state  and  the  cell 
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surface,  the  entire  simple  wave  can  be  approximated  by  a  single  char- 
acteristic.  In  the  same  way,  the  flow  on  the  cell  surface  must  be  given 
by  a  characteristic  coming  from  the  cell  on  the  other  side  of  the  cell 
surface;  thus,  the  fields  on  the  cell  surface  are  given  by  the  splution 
of  two  characteristics. 

i  ■ 

b.  Characteristic  Equations ' 

In  the  absence  of  nonuniformities  in  the  fields,  charac¬ 
teristic  properties  are  properties  whicfc  remain  constant  on  specific 
space-time  paths.  These  characteristic  properties  are  to  be  referred 
to  the  cell  surfaces;  consequently,  it  is  convenient  to  introduce  a  local 

rectangular  Cartesian  coordinate  system  on  each  cell  surface  as  depicted 
in  Figure  3. 


Figure  3.  Cell-Surface  Coordinate  System 


In  terms  of  this  coordinate  system,  it  is  shown  in 
Subsection  VM  of  Appendix  VI  that  the  hydrodynamic  characteristics 


rpa 


p  a 

a 


and  that  Ja±  is  constant  on  the  path 

n  =  n0+aa(Ma+i)t 

while  Ja_  is  constant  on  the  path 

n  =  nQ  +  a^Ma  -  l)  t 

whe- e  a,^  is  Jie  speed  of  sound  of  species  O  and 
Ma  = 


Similarly,  it  is  *hown  in  Subsection  VI-2  of  Appendix  VI 
that  the  radiation  transport  characteristics  are 

Jv±  ■  eRv±^QRvn  (49) 

and  that  Jv+  is  constant  on  the  path 

nQ+ct  fy/T  (50) 

while  Jv_  is  constant  on  the  path 

n0-ct//T  (51) 

where  QrVr  is  the  component  of  Qp^  along  m 

Maxwell's  equations  are  treated  in  Subsection  VI-3  of 
Appendix  VI.  It  is  shown  that  Maxwell's  equations  have  two  sets  of 
characteristics: 
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^2±  =  Bs2  *  ^S1 


Further,  Jj  and  j£  are  constant  along  the  path 


nQ  -  ct 


while  Ji_  and  J2+  are  constant  along  the  path 


n  +  ct 
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c*  Characteristic  Properties  on  Cell  Boundari 


The  characteristics  listed  in  Subsection  III-2b  define 
the  properties  on  the  cell  boundaries.  Thus,  let  cell  1  lie  on  the  nega¬ 
tive  n  side  of  the  cell  surface  and  cell  2  on  the  positive  n  side.  Then, 
if 


-1  <  Ma  <  1 


Equations  46  and  47  show  that 
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Equations  45  and  57  give 
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Assuming  tha„  the  integrand  in  Equation  59  is  nearly  constant  becaus< 
of  the  assumed  small  difference  between  properties  on  the  cell  surface 
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and  those  in  the  adjacent  cells. 
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Similarly,  Equations  45  and  58  lead  to 


<"an>, 


Fa 

,  “cell  surface 

cell  surface  *  (Paaah 


<ua">2  ‘  (6l) 


Solving  Equations  60  and  6i  for  the  cell  surface  properties. 
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The  density  on  the  cell  surface  Is  obtained  by  use  of 
Equation  397  of  Appendix  VI.  Thus,  if 


<«an>, 


cell  surface 


the  gas  on  the  coll  surface  is  deemed  to  have  come  from  cell  1,  and 
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Similarly,  if 

<Uan>cell  surface  <  0 
the  density  is  given  by 


^acell  surface  ~  ^0-2  + 


Pacell  surface  Pa2 
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With  Pa  and  Oa  known  on  the  cell  surface,  the 
remaining  thermodynamic  properties  on  the  cell  surface  are  found  by 
application  of  the  chosen  model  of  the  equation  of  state.  The  velocity 
tangential  to  the  cell  surface  is  taken  equal  to  the  tangential  velocity 
in  the  cell  from  which  the  gas  crosses  the  cell  surface. 

If  Equation  46  is  not  satisfied,  the  flow  is  super¬ 
sonic;  such  a  condition  means  that  the  gas  moves  faster  than  either 

characteristic.  In  this  case,  the  flowing  gas  drags  its  characteristics 
with  it.  Thus,  if 

Ma  ^  1  (68) 

the  flow  is  supersonic  from  cell  1  into  cell  2,  and  both  characteristics 
on  the  cell  surface  come  from  cell  1.  In  this  case,  the  flow  properties 
on  the  cell  surface  are  identical  to  those  in  cell  1 ; 
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Similarly,  if 
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the  flow  is  supersonic  from  cell  2  into  cell  1  and 
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If  the  cell  surface  itself  is  moving,  the  selection  of 
the  appropriate  set  of  cell  surface  properties  is  based  upon  the  motion 
of  the  characteristics  relative  to  the  moving  surface.  Since  the  char¬ 
acteristics  move  at  the  local  speed  cf  sound,  these  choices  may  be 
summarized  as  follows: 

I 

1.  If 


<“«„>,  -  -  wn>c 


th.  now  1,  'uperionic  from  cell  1  into  cell  2  and  Equations  69  through 
71  are  appi  ->priate; 

2.  If 

(77) 

the  flow  is  sv  per  sonic  from  cell  2  into  cell  1  and  Equations  73  through 
75  are  appropriate; 

3.  If  neither  of  inequalities  76  and  77  holds,  the  flow  is  subsonic 
relative  to  the  cell  surface,  and  Equations  62  through  67  are 
appropriate. 

The  characteristic  surface  properties  are  simpler  for 
the  equations  of  radiation  transport  because  these  characteristics 
propagate  at  a  speed  of  c//T;  consequently,  there  is  no  practical 
case  in  which  the  cell  surface  can  outrun  a  characteristic.  Thus, 
Equations  50  and  51  show  that 


^  v+^cell  surface 


(JvJ. 


^v-^cell  surface 


UvJ. 


(79) 


Inserting  Equation  49  into  Equations  78  and  79  and  solving  for 
the  cell  surface  properties  by  the  same  method  as  used  in  Equations 
60  through  63  , 
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Similarly  for  Maxwell's  equations.  Equations  54  and 


55  give 
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Using  Equations  52  , 

82  ,  and  83  and  solving  for  the 

surface  properties  gives 
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(86) 


(87) 


Similarly,  Equations  53  ,  84  ,  and  85  give 

,  (B‘Z>1+(B.2).t(E-l).-(E.,) 


^  S2^cell  surface 


(88) 
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3-  CONSTITUTIVE  EQUATIONS 

As  discussed  in  Subsection  Il-ld,  the  values  of  the 
collision  cross  sections  and  of  the  various  reaction  rates  will  be  assumed 
known.  The  values  initially  selected  will  be  obtained  from  a  cursory  check 
of  the  literature.  Since  the  purpose  of  the  initial  one- dimensional  calcu¬ 
lations  is  the  testing  of  the  sensitivity  of  the  numerical  solution  to  the 
various  parameters  of  the  problem,  such  a  choice  of  cross  sections  and 
reaction  rates  is  a  satisfactory  method  of  gaining  a  starting  point  for  a 
parametric  study.  A  more  thorough  evaluation  of  these  parameters  will 

have  to  be  undertaken  for  the  two-dimensional  overall- effects  code, 
however. 


The  detailed  study  of  initial  conditions  that  may  be 
required  for  predicting  the  debris  behavior  of  specific  bursts  is  being 
postponed,  at  least  until  the  two-dimensional  code  is  ready  for  running. 
If  the  parametric  studies  proposed  for  the  one- dimensional  code  indi¬ 
cate  a  lack  of  sensitivity  to  the  precise  initial  conditions,  then  such 
detailed  studies  may  not  be  warranted. 

a.  Source  Effects 

The  continuity  equation  for  species  a  requires  a 
knowledge  of  the,  volume  production  and/or  loss  rate  of  the  species. 

We  have  begun  to  investigate  the  various  cross  sections  and  reaction 
mechanisms  which  contribute  to  the  rates.  Among  these  are  ioniza¬ 
tion,  charge  transfer,  dissociation,  and  recombination  of  free  electrons 

with  positive  ions.  The  following  summarizes  the  preliminary  results 
of  our  survey  to  date. 


Capture  and  loss  of  electrons  by  fission  fragments 


Bell  (Reference  6)  and  Bohr  and  Lindhard  (Refer¬ 
ence  7)  have  shown  that  the  capture  and  loss  cross  sections,  ac  and 
al  •  respectively,  of  a  highly  ionized  particle  of  atomic  number  Z 
and  ionic  charge  z  moving  with  velocity  v  in  a  target  gas  of  moder¬ 
ately  high  atomic  number  Zt  are  given  by 


1/3  2 


3  2 


ac  “  Zt  z  (v0/v)  irao 

_  2/3  4/3  -3.  .  .2  2 

°l  ~  Zt  Z  z  (v/v0)  trao 

2 

Here  v0  is  c/137  and  irao  is  the  area  of  the  first  Bohr  orbit 
(8.  8x1 0”1 7  cm2). 


When  charge  equilibrium  is  attained  as  the  fission 
fragment  moves  through  the  gas,  the  capture  and  loss  cross  sections 
are  equal.  The  common  cross-section  under  equilibrium  conditions 
is  then  approximately 

/»»»»  vQ \  2 


•’  *  (ZZt)*/2(^)irao 

and  the  average  ionic  charge  is 
„  _  ^4/15  _  1/15/v ' 


(2)  Ionization  of  air  by  debri 


s  ions 


Perhaps  the  currently  most  useful  estimate  of  the 
cross  section  oi  for  ionization  of  a  gas  by  high-speed  ions  is  that  froi 
Firsov  (Reference  8).  His  result,  derived  for  atomic  systems,  is 

0i  =  -  ‘j  cm2  (94) 


3d 


where  v  is  the  relative  speed  of  the  colliding  pair  of  particles  and 
a0  and  v0  are  determined  by  the  relations 
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(95) 


ere,  Z,  and  Za  are  the  atomic  numbers  of  the  incident  and  target 
particles,  respectively,  and  Jo  is  the  smaller  of  the  ionization  energy 
(ev)  of  the  two  colliding  particles. 


Recent  measurements  by  Fite  et  al  (Reference  9) 
on  ionization  in  air,  jsfc  and  A-  gases  by  incident  beams  of  Al+  ions 
have  produced  the  data  shown  in  Figure  4.  Included  in  this  figure  are 
the  theoretical  cross  sections,  predicted  from  Equation  94  for  Al+ 
ions  in  argon  and  atomic  nitrogen.  It  is  seen  that  the  experimental 
curves  do  not  exhibit  quite  as  strong  a  velocity  dependence  as  does  the 
Fir sov  cross  section.  In  addition,  the  data  for  hfe  are  higher  than 


V  *  10  (cm/sec) 


Figure  4.  Ionization  Cross  Sections  as  a  Function  of  Velocity 
for  At  in  Ar ,  Ns ,  and  Air 


that  predicted  for  N.  It  seems  reasonable,  however,  that  in  dealing 
with  diatomic  molecules, the  value  of  o0  derived  for  atoms  should  be 
increased  somewhat.  This  would  then  bring  the  theoretical  and  exper 
mental  curves  more  into  line. 


(3)  Ionization  of  air  by  debris  electrons 

The  available  data  on  electron-induced  ionization  ir 
various  gases  have  recently  been  summarized  by  Kieffer  and  Dunn  (Ref¬ 
erence  10).  In  particular,  the  cross  sections  for  ionization  of  N,  O, 
I'fc.  Os,  NO,  N  and  Li+  by  electrons  with  energy  between  about  10  e' 
and  10  kev  are  available.  If  the  electrons  are  in  equilibrium  at  a  tem¬ 
perature  0e,  then  the  rate  coefficient  a  for  production  of  ions  M+ 
by  the  process 

e  +  M  -*  M+  +  2e  (96) 

is  given  by 


a  = 


<aiV> 


Here  is  the  electron  speed,  is  the  ionization  cross  section,  and 
the  brackets  indicate  an  average  over  all  speeds.  If  we  use  the  data 
presented  in  Reference  10,  together  with  an  assumed  velocity  distribu¬ 
tion  for  the  electrons,  Ot  as  a  function  of  electron  temperature  can  be 
obtained.  This  has  been  done  {Reference  11)  for  ionization  of  atomic 
oxYgen  and  nitrogen  by  assuming  a  Maxwell  velocity  distribution.  The 
results  are  shown  in  Figure  5.  Similar  results  can  easily  be  obtained 
for  the  other  species  by  numerical  integration  of  the  cross  section  data. 


Dissociative  ionization  of  air  mob 


Data  have  been  obtained  (Refer*-  je  10)  on  the  dis¬ 
sociative  ionization  proces s 
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Figure  5.  Ionization  Ratea  for  Atomic  Oxygen  and  Nitrogen 
Versus  Electron  Temperatures. 

e  +  Ms  -*  M  +  M+  +  2e  /OQ 

(Vo) 

for  electron  energies  between  about  10  ev  and  1  kev  and  product  ion 
(M  )  energies  greater  than  0.  25  ev.  Results  for  hfe  and  Cb  are 
shown  in  Figures  6  and  7. 

(5)  dissociation  of  ajr  molecules  bv  fa  afnm. 

Gerasimenko  and  Oksyuk  (Reference  12)  have  com¬ 
puted  the  cross  section  for  dissociation  of  diatomic  molecules  in  collisio 


log  10  E  (•«) 


Figure  6.  Cross  Sections  for  the  Dissociative  Ionization  of 

Molecular  Nitrogen  Yielding  Product  Ions  with  Kinetic 
Energies  Greater  than  0.  25  ev. 


Figure  7.  Cross  Sections  for  the  Dissociative  Ionization  of 

Molecular  Oxygen  Yielding  Product  Ions  with  Kinetic 
Energies  Greater  than  0.  25  ev. 


With  atoms.  It  turns  out  that  the  quantum  mechanical  expression  for 
the  cross  section  which  they  derive  gives  results  which  differ  at  most 
by  about  20  percent  from  the  classical  expression.  We  will,  there¬ 
fore,  adopt  the  simpler  classical  expression  for  the  dissociation  cross 
section. 


The  dissociation  cross  section  can  b 


e  written  as 


2  rT 

- 

j=l  •'D 


tfj(T)  w(T)  dT 


where  CTj(T)  is  the  cross  section  for  elastic  scattering  by  o-ie  of  the 
atoms  of  the  molecule,  T  is  the  energy  transferred  elastically  to  the 
atom,  D  is  the  dissociation  energy  and  w(T)  is  the  probability  that 
dissociation  will  occur.  Tm.  the  maximum  energy  transferred  elas¬ 
tically,  i8  given  in  terms  of  the  energy  E  of  the  incident  atom  and  the 


masses  M  and  Mj  of  the  incident  and  target  atoms,  by 


where 


4MMj 

^ _ J 

(M+M.)' 


(100) 


(101) 


The  summation  in  Equation  99  is  over  the  two  atoms  in  the  molecule. 

Classically,  w  =  1  ,  i.  e. ,  the  molecule  will  defi¬ 
nitely  dissociate  if  the  energy  transferred  elastically  to  one  of  the 
molecular  atoms  exceeds  the  dissociation  energy.  For  a  Coulomb 
interaction  potential,  we  have  (Reference  12), 


CTj(T)  dT 


2ire4Z2Zj24  dT 


M2  2 

Mj  v 


(102) 
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where  4  is  the  reduced  molecular  mass,  v  is  the  relative  speed 
between  incident  and  target  atoms,  and  Ze  and  Zje  are  respectively 
the  nuclear  charges  on  incident  and  target  atoms.  If  Tm  »  D,  as  will 
usually  be  the  case.  Equations  99  and  102  yield 
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It  is  not  difficult  to  use  the  more  accurate  screened- 
Coulomb  interaction  or  the  Thomas -Fermi  function  to  obtain  a  slightly 
more  accurate  expression  for  the  cross  section. 


In  addition  to  the  collisional  effects  discussed  above, 
the  decay  rates  from  thermodynamic  nonequilibrium  to  equilibrium  are 
expected  to  be  important  at  higher  altitudes.  These  decay  equations  are 
generally  given  in  the  form 


dtp 

dt  “  wtyi  »  $ 3  <  •  •  • »  T)  (104) 


where  cpx ,  Cfc ,  ...,  cpN  are  the  N  properties  which  are  out  of  equilib- 
rium  and  T  is  t^ie  temperature.  The  function  U)  in  Equation  104 
frequently  is  uighly  nonlinear,  and  recourse  must  be  had  to  numerical 
means  of  integrating  the  equation.  These  rate  equations  frequently  have 
very  large  time  constants,  however,  and  conventional  Runge-Kutta  tech¬ 
niques  therefore  do  not  work  well.  However,  a  modification  of  a  method 
suggested  by  Certaine  (Reference  13)  has  been  developed;  this  modifica¬ 
tion  appears  accurate  and  efficient  in  terms  of  computer  time. 


(!) 


a 


The  function  to  in  Equation  104  is  approximated  by 


-Dcpc  +  @(t) 


(105) 
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where 


Equation  111  is  computed  in  three 


(1)  The  approximations 


D 


and 


(112) 


(113) 


are  used,  and  a  preliminary  value  of  cpa(At)  is  computed; 

(2)  Using  the  preliminary  value  of  ®  (At),  final  values  of  D  and 
cPoj^At)  are  computed; 

(3)  The  final  values  of  D  and  %^t)  are  u.ed  in  Equati<)n  (o 
obtain  the  final  value  of  cpa(At). 


Equation  111  has  been  coded  and  tested  by  being 
employed  to  calculate  the  relaxation  to  equilibrium  of  a  mixture  of 
Oz,  Qs  ,  O,  and  O  with  initial  conditions  of  10,  000  °K  and 

100  percent  0+.  The  computed  relaxation  was  checked  against  results 

obtained  by  Runge-Kutta  techniques  and  found  satisfactory. 
c<  Correctness  of  Treatment 

The  treatment  of  the  constitutive  equations  as  dis¬ 
cussed  in  Subsection  III- 3a  is  decoupled  from  the  computation  of  the 
fluxes  crossing  the  cell  surfaces  as  discussed  in  Subsection  III-2. 

The  correctness  of  this  decoupling  is  seen  by  recalling  that  the  sur¬ 
face  fluxes  have  to  do  with  the  transport  of  mass,  momentum,  and 
energy  across  the  cell  surfaces.  As  pointed  out  by  Chapman  and 
Cowling  (Reference  1),  however,  the  transport  of  molecular 
properties  in  a  rarefied  gas  is  caused  almo  entirely  by  the  free  motion 
of  the  particles  between  collisions  and  only  t .  ligibly  to  the  transf 
at  collisions  over  the  distance  separating  the  two  colliding  particles. 
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Thus,  for  the  rarefied  gases  oon.idered  here,  partiole  collisions  and 
molecular  relaxations  can  have  no  effect  upon  the  surface  fluxes. 

Consequently,  the  decoupling  of  the  collisional  effects  from  the  sur- 
face  flux  calculations  is  valid. 

4‘  BOUNDARY  CONDITIONS 

The  appropriate  types  of  boundary  conditions  and 
the  means  by  which  they  may  be  computed  follow  directly  from  the 
cell  surface  properties  given  in  Subsection  III-2c.  The  boundary  con¬ 
ditions  must  he  such  that  these  surface  properties  can  be  computed; 

generally,  this  requirement  means  that  all  of  the  fields  outside  of 
the  boundary  must  be  prescribed. 


In  the  case  of  the  hydrodynamic  equations,  how¬ 
ever.  an  additional  type  of  boundary  condition  arises  if  the  flow  is 

supersonic  out  across  the  boundary.  In  this  case.  Equations  69 

through  71  show  that  the  surface  properties  on  the  cell  boundary 
ere  given  by  tne  properties  in  the  last  cell  within  the  boundary. 
Thus,  the  external  fields  need  not  be  prescribed  for  this  case. 


n  a  dition,  it  is  possible  to  replace  the  surface 
property  calculations  of  Subsection  111.2c  by  specified  surface  proper, 
ties.  For  example,  if  a  cell  surface  is  known  to  coincide  with  a 
hydrodynamic  shock  wave,  a  subroutine  containing  the  Rankine- 

Hugoniot  relations  can  be  used  for  the  calculation  of  the  cell  surface 
properties. 


5.  _MESH  GEOMETRY 

The  finite  difference  equations  of  Subsection  III-l 
have  been  written  in  a  generalised  geometry  which  facilitates  the 
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aligning  of  the  cell  sides  with  any  discontinuities  in  the  flow.  This 
generalized  cell  geometry  coupled  with  the  motion  of  the  mesh  implies 

that  the  cross  sections  of  the  cells  are  quadrilaterals,  but  not  necess¬ 
arily  rectangles. 

a*  One  or  Two  Dimensions 

Three  coordinate  systems  are  used;  plane,  cylin¬ 
drical,  and  spherical.  In  the  first  two  of  these,  computations  may  be 
chosen  to  be  either  one-  or  two-dimensional;  calculations  in  the  spheri¬ 
cal  system  must  be  one -dimensional.  The  ability  to  choose  between 
one-  and  two-dimensional  calculations  results  from  the  surface  flux 
technique  of  differencing  the  equations;  a  one- dimensional  calculation 
is  achieved  if  two  of  the  cell  sides  are  dropped  from  the  summations 
in  Equationr  28  through  30  ,  36  through  38  ,  and  41  through  44  . 
This  dropping  is  achieved  by  setting  the  flag  N DIMEN  equal  to  unity. 

In  addition,  the  z- component  equations  may  be  neglected  in  some  of 
these  calculations  (but  need  not  be).  The  resulting  calculations  will 
be  executed  as  quickly  as  if  the  code  had  been  written  specifically  for 
one -dimensional  calculations. 

b.  General  Quadrilateral  Shapes 

The  general  quadrilateral  cross  section  of  a  cell 
is  depicted  in  Figure  2.  In  the  case  of  cylindrical  geometry,  this 
cross-sectional  area  is  revolved  around  the  z  axis  (Figure  1)  to  gen¬ 
erate  the  computational  cell.  In  the  case  of  spherical  geometry,  the 
cell  is  a  spherical  shell,  as  illustrated  in  Figure  1.  The  choice 
between  these  three  geometries  is  made  by  setting  the  geometry 
option  word,  IGEOM,  equal  to  one  for  plane  geometry,  two  for  cylin¬ 
drical  geometry,  and  three  for  spherical  geometry. 
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c*  Method  of  Specification 

The  parameters,  needed  to  specify  the  cells  are 
epicted  in  Figure  2.  Each  cell  is  bounded  by  two  rays  and  two 
arcs;  cell  {i+T/2.  jd-1/2,  is  bounded  by  rays  i  and  i+1,  and  by 

arcs  j  and  j+1  .  These  rays  and  arcs  intersect  to  form  the  four 

corner  points  of  the  cell:  (i.  j),  (i+1,  J),  (i,  j+1),  and  (i+1,  j+1). 

The  totality  of  cells  between  any  two  rays  is  referred  to  as  a 
column  of  cells. 

The  arcs  constitute  the  moving  portion  of  the  mesh. 

As  a  result  of  this  motion,  the  arcs  may  become  discontinuous  as 
shown  in  Figure  8. 
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Figure  8.  Discontinuous  Arcs 

The  rays,  however,  are  fixed  straight  lines  during  any  one  calcu¬ 
lation.  Thus,  cell  (i+l/2,j+l/2)  is  completely  specified  if  the 
ngles  a.  and  a.+  1  of  the  adjacent  rays  are  given,  if  the  intercept! 

zi  and  zitl  of  the  rays  on  the  z  axis  are  given,  and  if  the  dis- 

tances  d.  . ,  d  d  j 

i.  j+1  *  i+1,  j  and  di+l,j+l  along  rays  i  and  i+1 

from  the  z  axis  to  the  corners  of  the  cell  are  given. 
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As  can  be  seen  from  Figure  2,  all  of  the  areas  and 
volumes  associated  with  a  cell  are  computable  in  terms  of  triangles 
and  triangles  of  revolution  about  the  z  axis.  Thus,  the  development 
of  the  formulae  /or  the  cell  geometry  is  straightforward.  The  results 
are  summarized  in  Appendix  VII. 

The  first  ray  (i=l)  and  the  last  ray  (i=l)  are  consid¬ 
ered  to  be  external  boundaries,  and  have  their  surface  fluxes  pre¬ 
scribed  in  subroutines.  In  this  fashion,  the  changing  of  boundary 
conditions  is  reduced  to  the  changing  of  a  few  subroutines. 

In  the  same  way,  the  first  arc  (j=l)  and  the  last  arc 
^=Ji+l/2*  are  considered  to  be  external  boundaries.  In  addition,  the 
moving  mesh  can  accommodate  itself  to  moving  boundaries;  thus,  any 
number  of  internal  arcs  may  be  aligned  with  discontinuities  in  the 
internal  flow.  For  example,  a  contact  surface  might  be  represented 
by  arc  j=5  in  column  (i+1/2),  arc  j=6  in  column  (i+3/2),  and  so 
forth.  All  of  these  arc  boundaries,  both  external  and  internal,  have 
their  surface  fluxes  prescribed  in  subroutines;  thus,  external  condi¬ 
tions  and  internal  discontinuities  are  easily  changed. 

In  any  given  column,  the  arcs  which  are  not  associated 
with  either  external  boundaries  or  internal  discontinuities  are  positioned 
along  the  rays  in  such  a  way  as  to  provide  equal  spacing  between  the 
discontinuities.  Thus,  the  location  of  the  points  of  intersection  of  all 
non-boundary  arcs  with  the  rays  can  be  computed  if  the  points  of  inter¬ 
section  of  the  boundary  arcs  are  given.  For  example,  if  arc  j=J.+  1/2 

has  intersections  with  rays  i  and  i+1  located  at  d  and  d 

....  i*  J  i+ 1 ,  J  ’ 

and  it  the  next  boundary  arc  is  an  internal  boundary  on  arc  j=N  with 

points  of  intersection  at  d.^  N  and  d.+  ^  N;  then  the  intermediate 

points  of  intersection  are  located  at 
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d  +  1>  H  J 
i»j+l  J  -  N 


(114) 


i+ 1 »  j 


d  +  -  i+l«  N  di+l,  J 

i+ 1  j  j+ 1  J  -  N 


(115) 


d.  Multiple  Boundaries 

As  can  be  seen  by  an  examination  of  Figure  8,  the 

moving  mesh  may  cause  the  face  on  ray  i  of  cell  (i+1  /2,  j+1/2)  to 

share  ray  areas  with  the  face,  of  several  cells  in  the  next  column. 

In  such  a  case,  the  characteristic  property  calculations  of  Subsection 
HI- <ic  are  carried  out  on  each  portion  of  the  face,  and  the  summations 
in  the  integrated  equations  of  Subsection  III-l  are  taken  over  all  portions 
of  the  cell  face.  The  ray  area  formulae  in  Appendix  VII  can  be  used 
for  the  computation  of  each  of  these  partial  cell-face  areas. 

This  technique  of  summing  the  fluxes  on  multiple  bound¬ 
aries  has  been  used  successfully  in  References  3  and  5.  The  method 
ensures  the  conservation  of  mass,  momentum,  and  energy. 

(>•  MESH  MOTION 

To  follow  the  motion  of  the  external  and  internal 
boundary  arcs  (if  they  move),  it  is  necessary  to  allow  all  of  the  arcs  to 
move.  Thus,  the  motion  of  the  boundary  arcs  is  determined  by  the  kine¬ 
matics  of  the  boundary,  and  the  motion  of  the  remaining  arcs  is  scaled 
to  maintain  equal  arc  spacing  between  boundaries. 

a.  Node  Velocities 

The  nodes  are  the  points  of  intersection  between  the  arcs 
and  the  rays.  Those  nodes  which  lie  on  arcs  corresponding  to  internal 
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or  external  boundaries  have  their  velocities  determined  by  the  velocity 
of  the  arc  normal  to  itself.  The  means  of  computing  these  node  veloc¬ 
ities  can  be  explained  with  the  aid  of  the  geometry  depicted  in  Figure  9: 

r 


The  geometry  illustrated  in  Figure  9  is  easily  calculated  from  the  for¬ 
mulae  of  Appendix  VII.  The  arc  velocities.  Wl  and  Wa ,  are  given 
by  the  kinematic,  of  the  boundary;  for  example,  a  contact  surface  moves 
at  the  local  material  velocity  and  a  shock  wave  moves  a,  that  velocity 
winch  renders  the  pressure,  density,  and  velocity  jumps  across  the  dis- 
contmuity  compatible  with  the  Rankine-Hugoniot  equations.  These  kine- 
matrc  velocity  calculations  are  included  in  the  internal  or  external  boundary 
subroutines  appropriate  to  the  arc  in  question. 


In  terms  of  the  definitions  illustrated  in  Figure  9,  the  node 

velocity  is 

v_  =  QiWg / sin  Ys  +  -ta Wi / sin  Yl"| 

L  ir+~u  j  (H6) 

The  velocities  of  the  nonboundary  nodes  along  a  ray  are 
scaled  to  maintain  equal  node  spacing,  between  boundary  nodes.  For 

this  purpose.  Equations  114  and  115  are  used  with  the  distances 
replaced  by  the  node  velocities. 
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b. 


Cell  Surfaces 


The  surface  faces  of  the  C3lls  on  the  fixed  rays  do  not 
move,  but  the  faces  on  the  arcs  do  move.  The  velocity  of  an  arc  face 
normal  to  itself  is  evaluated  by  calculating  the  volume  swept  out  by 
the  arc  face  during  the  time  step  and  dividing  this  volume  by  the  product 
of  the  face  area  and  the  time  step.  This  procedure  is  used  to  ensure 
that  volume  is  conserved  by  the  floating  mesh.  Thus,  if  the  position  of 

node  (i,  j)  at  the  start  of  the  time  step  is  d.^  (0) ,  the  position  at  the 
end  of  the  time  step  is  * 


d-  -(0)  +  v.  .  •  At 
J  J 


in  n 


where  v_  is  the  velocity  of  the  node  along  ray  i  as  given  by  Equa¬ 
tion  114  or  Equation  116.  The  node  positions  before  and  after  the 

time  step  are  used  in  the  cell  volume  formulae  of  Appendix  VII  to 
obtain  the  volume  swept  out  by  the  cell  face  during  the  time  step,  AV. 
The  velocity  of  face  (i+1/2,  j)  is  then  found  from 


Wj  =  -Av/(wr*> 


(118) 


In  previous  calculations  with  this  type  of  moving  mesh 
(Reference  5),  it  has  been  found  necessary  to  recompute  the  velocity 
of  any  internal  boundary  arcs  according  to  Equation  118  for  the  sur¬ 
face  flux  calculations;  the  velocity  so  computed  differs  only  slightly 
from  the  velocity  given  by  the  kinematic  subroutine  corresponding  to 
the  arc,  but  this  slight  difference  is  sufficient  to  cause  inaccuracies 
because  of  lack  of  conservation  of  volume. 

7-  STABILITY  ANALYSIS  AND  TIME-STEP  CALCULATION 

The  stability  analyses  used  on  these  equations  is  based 
upon  the  von  Neumann  necessary  condition  as  presented  by  Ric.itmyer 
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(Reference  14).  The  equations  are  linearized,  the  amplification 
matrix  is  found,  and  the  eigenvalues  of  the  amplification  matrix  are 
determined.  Stability  requires  that  these  eigenvalues  be  less  than  unity. 


Two-Dimensional  Stability 


It  is  shown  in  Subsection  3.  J  of  Appendix  3  of  Reference 
3  that  the  two-dimensional  code  will  be  stable  if  the  time  step  is  chosen 
to  satisfy  the  inequality 


At  At 

<  — L_ 

“  Atr  +  Atz 


(119) 


where  Atr  and  Atz  are  the  time  steps  for  which  the  one- dimensional 
calculations  in  the  radial  and  axial  directions,  respectively,  are  stable. 


One- Dimensional  Stabiliti 


For  the  hydrodynamic  equations,  it  is  shown  in  Appendix 
G  of  Reference  3  that  the  numerical  method  given  here  is  stable  in 
the  radial  direction  if  the  time-step  Atr  satisfies  the  inequality 


(120) 


in  every  cell  of  the  mesh  and  for  each  species.  Here,  hr  is  the  dis¬ 
tance  across  the  cell  in  the  radial  direction,  ur  is  the  velocity  in  the 
radial  direction,  and  a  is  the  speed  of  sound  within  the  cell.  A  simi¬ 
lar  expression  holds  for  the  axial  time  step,  Atz.  Equation  119  is 
valid  if  the  cell  surface  properties  of  Subsection  3.  2.  3  are  evaluated 
at  the  start  of  the  time  step.  It  is  found  in  Reference  3  that  the  time 
step  permitted  by  Equations  119  and  120  is  sufficiently  long  to  allow 
large  percentage  changes  (as  much  as  50  percent)  of  the  flow  properties 


within  a  cell.  Since  larger  percentage  change,  would  raise  the  ques¬ 
tion  of  accuracy,  this  time-step  limit  is  considered  satisfactory. 

For  Maxwell's  equations  and  the  equations  of  radiation 
transport,  however,  the  situation  is  somewhat  different.  If  the  cell 
surface  properties  are  evaluated  at  the  start  of  the  time  step,  it  can 
be  shown  that  the  stable  time  step  must  satisfy  the  inequality 


4*r  <  hr/c 


(121) 


where  c  is  the  speed  of  light.  Obviously.  Equation  121  presents 
far  too  stringent  a  limitation  on  the  time  step  for  a  practical  numerical 
procedure.  However,  if  the  cell  surface  properties  are  evaluated  at 
the  end  of  the  time  step,  it  can  be  shown  (Appendix  VIII)  that  the  numer¬ 
ical  procedure  is  stable  for  any  time  step.  Consequently,  this  option  of 
evaluating  .he  surface  properties  at  the  end  of  the  time  step  is  chosen 
for  Maxwell's  equation,  and  the  equations  of  radiation  transport.  In 
the  case  of  one-dimensional  calculations,  the  resulting  numerical  equa¬ 
tions  are  easily  inverted  and  solved  directly,  but  in  the  case  of  two- 
dimensional  calculations,  an  implicit  set  of  equations  results.  Although 
there  exists  a  number  of  relaxation  methods  for  solving  implicit  equa- 
tions,  no  method  has  been  selected  yet. 

8*  order  OF  COMPUTATIONS  IN  THE  CODE 

The  general  flow  of  logic  in  the  code  is  indicated  in  detail 
in  the  flow  charts  in  Appendix  II.  The  various  formulae  presented  in  this 
section  are  coded  into  subroutines  and  called  as  needed.  The  flow  charts 
not  only  show  the  flow  logic,  but  also  give  the  FORTRAN  II  coding  used. 

The  code  does  not  yet  contain  Maxwell's  equations  or  the 
equations  of  radiation  transport.  Furthermore,  the  code  has  not  yet 
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bCen  C°mpletely  checked  <>"*.  although  some  test  calculations  have 
been  run  successfully. 


charts. 


A  few  minor  subroutines 


are  not  shown  in  the  flow 


SECTION  IV 
DISCUSSION 


A  numerical  technique  has  been  developed,  but  not 

yet  fully  tested,  for  the  study  of  the  dynamics  of  high  altitude  nuclear 
bursts. 

1.  METHOD  OF  USE 

The  numerical  technique  developed  here  can  be  used 
for  both  one -dimensional  and  two-dimensional  studies. 

a*  One -Dimensional  Coupling  Studies 

It  ts  proposed  to  use  the  physics  code  described  above 
in  a  series  of  one-dimensional  calculation,  to  study  the  coupling  to  the 
atmosphere  of  the  debris  from  a  high  altitude  explosion.  Separate 
studies  are  proposed  for  colUsional  coupling,  electromagnetic  pickup, 
and  pickup  by  scattering  off  of  magnetic  turbulence. 

The  collisional  coupling  studies  will  be  carried  out  by 
use  of  expressions  for  the  rate  of  change  of  momentum  caused  from 
collisions.  These  expressions  are  functions  of  the  collision  cross  sec¬ 
tions;  it  is  proposed  to  vary  the  densities,  energy  of  the  debris,  and 
cross  sections  within  their  range  of  uncertainty  to  observe  the  transi¬ 
tion  from  collisionless  flow  to  fully  coupled  flow.  Scaling  laws  for  this 
transitional  region  can  then  be  developed. 

The  two- stream  instability  also  can  be  studied  with  the 
one-dimensional  physic,  code.  Two  neutral  streams  of  ions  and  electrons 


are  allowed  to  approach  one  another  in  the  absence  of  colli sional 
terms.  The  presence  of  the  complete  Maxwell's  equations  in  the  code 
ensures  that  any  electric  field  generated  by  density  perturbations  will 
be  computed.  This  type  of  computation  should  give  information  as  to 
the  effective  role  of  the  two-stream  instability  in  the  coupling  process. 
Here  again,  scaling  laws  giving  the  magnitude  of  the  perturbation  and 
the  degree  of  pickup  can  be  developed. 

The  one-dimensional  code  also  can  be  used  for  studies 
>1  o  Tiling  by  scattering  jf  ions  off  of  magnetic  irregularities.  The 
magnetic  field  is  assumed  to  be  trapped  in  the  hydrodynamic  turbulence 
in  accordance  with  the  model  of  Reference  15.  The  turbulence  itself 
is  computed  by  introducing  a  turbulent  viscosity  coefficient  into  the  code, 
and  using  the  scaling  laws  of  Kolmogorov's  theory  (Section  32,  Refer¬ 
ence  16).  The  cross-section  for  scattering  is  given  in  Reference  15. 

By  coupling  the  scattering  off  of  the  magnetic  irregularities  to  the 
general  flow  by  means  of  the  physics  code,  it  will  be  possible  to  de¬ 
duce  whether  or  not  such  a  phenomenon  can  cause  coupling. 

b.  Two-Dimensional  Studies 

The  physics  code  discussed  above  also  has  been  developed 
in  a  two-dimensional  form  (both  plane  and  cylindrical  symmetries).  It 
is  proposed  to  use  this  form  of  the  code  to  study  "Longmire"  coupling 
and  to  compute  the  gross  structure  of  the  expanding  debris. 

The  concept  of  "Longmire"  coupling  is  fundamentally 
two-dimensional,  because  it  involves  the  generation  of  a  magnetic  field 
by  a  stream  of  ions  and  a  turning  of  additional  ions  by  the  magnetic 
field.  The  ability  to  remove  the  collisional  terms  from  the  code  plus 
the  presence  of  Maxwell's  equations  in  the  code  permits  the  simulation 
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of  this  coupling.  The  computation  is  set  up  in  cylindrical  coordinates. 
A  small  region  of  ionized  particles  in  the  center  of  the  computational 


mesh  is  given  an  initial  outward  velocity,  and  the  ambient  conditions 
outside  this  expanding  region  are  taken  as  those  of  an  ionized  gas  at 
rest  in  a  magnetic  field.  The  initial  expansion  energy,  degree  of  ioniza¬ 


tion,  and  initial  density  can  be  varied  in  a  parametric  study  of  this  form 

of  coupling.  If  coupling  does  occur,  appropriate  scaling  laws  can  be 
deduced. 


sc 


aung  laws  xor  the  various  type 


-  J  L. - —  lldVC 

been  deduced,  it  will  be  possible  to  use  the  two-dimensional  code  for 
studies  of  the  overall  expansion  of  the  debris  from  high  altitude  explosions. 
These  scaling  laws  can  be  introduced  at  the  expanding  outer  boundary  of 
the  debris  in  the  same  fashion  as  the  Rankine-Hugoniot  equations  are 
introduced  at  the  outer  boundary  of  conventional  hydrodynamics  codes. 


It  is  worth  noting  that  the  numerical  procedure  developed 
here  can  be  used  either  as  an  Eulerian  procedure  or  as  a  Lagrangian 
procedure.  In  a  one-dimensional  calculation,  for  example,  the  code 
will  be  Eulerian  if  the  mesh  is  fixed  in  space;  on  the  other  hand,  the 

code  will  be  Lagrangian  if  each  cell  surface  is  caused  to  move  at  the 
local  fluid  velocity. 


2-  FURTHER  development 

The  proposed  physics  code  is  capable  of  yielding  detailed 
numerical  information  on  the  problem  of  the  coupling  of  the  debris  of 
high  altitude  nuclear  explosions  to  the  ambient  atmosphere.  The  code 
is  also  capable  of  incorporating  this  information  in  a  manner  suitable 
for  the  computation  of  the  late-time  debris  structure.  Most  of  the 
code  as  described  in  the  loregoing  section  has  been  written  but  not 
debugged.  Maxwell's  equations  have  not  yet  been  coupled  into  the  code, 


59 


s  equations  in  two 


and  the  theory  for  the  treatment  of  Maxwell' 
dimensions  is  not  yet  complete. 

For  the  one -dimensional  code,  therefore,  it  is  pro¬ 
posed  that  the  first  task  be  the  debugging  of  the  existing  code.  There 
are  many  solutions  of  channel  flows  which  may  be  used  to  demonstrate 
the  accuracy  of  the  code.  The  second  proposed  task  is  the  development 
of  subroutines  for  the  various  thermodynamic  and  kinetic  models.  It 
is  suggested  that  an  equilibrium  air  equation  of  state  be  coded  into  a 
subroutine,  and  that  a  set  of  reaction  rate  equations  for  the  decay  from 
a  nonequilibrium  state  to  the  equilibrium  state  be  developed  and  coded. 

In  addition,  ic  is  suggested  that  this  second  task  include  the  development 
of  subroutines  which  give  the  source  of  momentum  due  to  particle  col¬ 
lisions  and  due  to  scattering  off  of  magnetic  irregularities.  These  last 
two  subtasks  will  involve  finding  or  developing  mathematical  expressions 
for  the  cross-sections  for  the  proposed  reactions. 

The  third  task  proposed  for  the  one-dimensional  code 
is  the  adding  of  Maxwell's  equations  to  the  code  and  the  subsequent  de¬ 
bugging.  The  theory  of  the  characteristic  flux  method  as  applied  to 
Maxwell' s  equations  in  one  dimension  is  completely  developed,  and  the 
relevant  stability  analyses  have  been  made.  These  equations  can  be 
coded  in  a  fashion  which  is  stable  for  any  time  step. 

For  the  two-dimensional  code,  it  is  proposed  that  the 
first  task  be  the  debugging  of  the  existing  code.  The  blast  penetrat:  >n 
calculations  of  Reference  3  can  be  used  as  one  standard  during  the  de¬ 
bugging,  and  the  point  explosion  results  given  in  Sections  11  and  12 
of  Chapter  IV  of  Reference  17  can  be  used  as  a  second  standard.  These 
proposed  standards  cover  the  cases  of  strong  interactions  of  high  energy 


streams  with  associated  material  interfaces  and  of  high  energy 
release  in  a  small  region. 


The  second  task  proposed  for  the  two-dimensional 
code  is  the  development  and  introduction  of  the  appropriate  thermo¬ 
dynamic  routines.  These  routines  would  include  all  of  those  developed 
for  the  one-dimensional  code,  plus  all  of  the  scaling  laws  for  coupling 
developed  by  use  of  the  one-dimensional  code.  The  latter  class  of 


thermodynamic  routines  would  serve  as  proven  coupling  law.  which  can 
be  used  in  the  computation  of  the  overall  structure  of  the  expanding  debris 


Finally,  it  i.i  proposed  that  a  method  be  developed  for 
coupling  Maxwell's  equations  into  the  two-dimensional  code.  The  charac¬ 
teristic  flux  technique  has  been  applied  to  these  equations,  but  the  resulting 
stability  matrix  is  not  explicitly  solvable  and  has  not  been  evaluated  numer¬ 
ically.  It  is  suggested  that  evaluation  of  the  stability  matrix  be  undertaken, 
and  that  the  equations  be  added  to  the  code  if  proven  stable.  If  the  equation’. 


are  not  stable,  it  is  suggested  that  attempts  to  generate  a  stable  technique 
be  undertaken. 
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APPENDIX  I 


CONSERVATION  EQUATIONS  FOR  A  REACTING 
NONEQUILIBRIUM  MIXTURE  OF  GASES 


The  derivation  of  the  usual  hydrodynamic  conservation 
equations  for  a  ncnreacting  mixture  of  gases  from  Boltzmann's  equa¬ 
tions  is  given  in  many  texts  on  physical  gas  dynamics.  A  similar 
derivation  is  presented  here  in  order  to  demonstrate  and  document 

the  consistence  cf  the  definitions  used  in  Sections  II  and  III  of  this 
report. 


1. 


DEFINITION  OF  KINETIC  PARAMETERS 


Assume  a  chemically-reacting  mixture  of  gases,  and 
assume  that  all  interactions  are  weak;  i.  e. .  no  electromagnetic  or 
radiative  collision  effects.  Assume  N  species,  including  electrons 
(the  species  may  be  ionized). 


Let  pa  denote  the  mass  per  unit  volume  of  species  a 
and  let  p  denote  the  mass  of  mixture  per  unit  volume: 

N 

P  SP<X  (122) 

Let  fa(u,  r,  t)  be  the  single-particle  distribution  function  for  species 
a  and  let 

na  =  number  density  of  particles  of  species  a 

ma  “  mass  of  a  particle  of  special  a 


Then 


'  »■ 

1 


na(r,  t) 

=  ff  du 

*  a  a 

(123) 

and 

=  mana(7,  t) 

(124) 

as 

The  mean  velocity  of  particles  of 

species  a  is  defined 

<%> 

fuf  du 

J  CC  CC  CL 

f{ ctd% 

=  —  fu  f  du 

n  »  a  a  a 
a 

(125) 

The  mass  velocity  of  the  mixture. 

u,  is  defined  as 

u 

_  y'  . 

a=l 

(126) 

is  defined  as 

The  thermal  velocity  of  a  particle  of  species  a,  va. 

-* 

V 

a 

*  V  <v 

(127) 

Note  that  Equation  127  is  not  the  conventional  definition  of  thermal 
velocity,  which  is  expressed  in  terms  of  the  mass  velocity  of  the 
mixture.  It  is  felt  that  for  the  situations  considered  here,  which  are 
nowhere  near  being  equilibrium  mixtures,  species  temperatures  are 
better  based  on  the  motion  of  the  individual  species  than  on  the  motion 
of  the  mixture.  The  latter  motion  is  given  by  the  mean  thermal  vel¬ 
ocity  of  a  particle  of  species  a,  <va>,  defined  as 

=  ua-u  (128) 
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obviously. 


;a-<v 


5-<%> 


(129) 


The  diffusion  velocity  of  particles  of  species  a,  V,  is  usually  defined 


<  V  -  u 


(130) 


Thus,  Equation  129  becomes 


v  -  <  v  > 

a-  a 


(131) 


2‘  BOLTZMANN  EQUATION  FOR  SPECIES  r. 


Using  tensor  notation,  the  time- rate -of- change  of  specie* 
a  within  a  fixed  volume  in  (u,  r)  spaqe  is 


^(»oV) 


*  3  ffc  ad"d? 

=  d3dr 


(132) 


This  change  results  from  two  physical  phenomena:  i)  particle  flux 
across  the  surface  of  V;  and  2)  particle  sources  within  V.  The 
change  arising  from  particle  fluxes  is 


*(naV»flux 


-/fct'*’1'V 

e  «*• 


(133) 


where 


generalized  velocity  of  a  particles  in 
(%»  ra)  space, 

unit  outward  normal  to  S, 


03 


and 
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(139) 


Using  Equations  136  and  137  in  Equation  134  gives 

3XV)fluX  *  idV 

a.  Source  Terms 

The  particle  sources  within  the  six- dimensional  volume 

v  are  (1)  creation  of  particlee  within  physical  space  by  chemical  or 

nuclear  processes;  and  (2)  creation  of  particles  within  velocity  space 

by  scattering  from  another  element  of  velocity  space.  Note  that  the 

nuclear  processes  involved  here  may  emit  or  absorb  photons,  and 

the  numerical  processes  being  developed  for  these  equations  can  treat 

photons  as  well  as  other  species.  Let  the  particle  sources  be  defined 
as  follows: 

Ra  =  ^ber  of  particles  of  species  a  generated 

per  unit  time  in  unit  phase  space  volume  around 
point  (u,  r)  by  all  processes; 

Ra  =  number  of  particles  of  species  a  lost  per  unit 
time  in  unit  phase  space  volume  around  point 

/-4 

(u,  r)  by  all  processes. 

Then 

dt  ^sources  =  dY  (140) 

But 


sources 


(141) 


Combining  Equations  132  ,  139,  140  and  141, 


(142) 
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Noting  that  the  volume  of  integration  in  Equation  142  is  arbitrary, 
the  Boltzmann  equation  results: 


3f 

a 

dt 


dr 


Rq 


R 


a 


(143) 


The  source  rates  in  Equation  143  are  net  rates  result¬ 
ing  from  all  possible  reactions;  thus. 


R 

a 

y '  <2)  y  <3) 

tr <144) 

where 

(2) 

rate  at  which  particles  of  species  p  and  y 

react  in  two-body  collisions  to  form  particles 

(3) 

of  species  a,  and 

Ra,  3y;6 

rate  at  which  particles  of  species  (3  and  y 

react  in  three-body  collisions  (in  the  presence 
of  species  6)  to  form  particles  of  species  a 

Equation  144  obviously  can  be  generalized  to  include  higher  than 

three-body  collisions,  but  such  effects  seem  unimportant. 


There  is  a  similar  equation  for  R  : 

v-<3> 

rate  at  which  particles  of  species  a  and  (3 
react  in  two-body  collisions  to  annihilate  par¬ 
ticles  of  species  a,  and 


R 


a 


where 


(2) 


R 


a  p 
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and 


_(3) 

Rag,  y  "  rate  at  which  particles  of  species  a  and  g  react 
in  three-body  collisions  (in  the  presence  of  species 
Y)  to  annihilate  particles  of  species  a. 

It  should  be  remembered  that  the  terms  -form"  and  -annihilate"  as 
used  in  the  above  definitions  refer  to  an  element  of  six- dimensional 

phase  space,  and  include  scattering  from  one  element  of  velocity  space 
to  another  element  of  velocity  space. 


The  above  source  rates  ca,  be  expressed  in  terms  of 

cross  sections  as  follows.  Suppose  that  a  particle  of  species  a  and 

a  particle  of  species  8  traveling  with  relative  velocity 


u  -  u 
1  a  3 


(146) 


have  a  cross  section  for  the  reaction 


a»  P  “*  Y»  6,  e 


(147) 


nfe  ,  Ti  ,  VI  ;  g  ) 
'  Y  6*  s  8a g/ 


(148) 


Where  the  cross  section  depends  upon  the  velocities  of  the  product 

particles.  Then  the  rate  of  formation  of  particles  y.S.e  because  of  this 

particular  two-body  interaction.  T  is 

yf>e,  ag  8 

rY6e,  ag(uY*  V  V  t)  =  ?.  t)  •  fp(up,  r,  tj 

'  gag’  4y  V  V*fcg)dSadSg  (149) 

The  rate  of  formation  of  any  one  species  by  two-body  collisions  is 


(2)  , 

RY,agVr'V  =  .//rY6e,  aivVV'*1)^^  050) 


b9 


It  is  clear  that  if  a  standard  form  (such  as  the  Maxwellian)  for  the 
distribution  functions  f^  and  f^  can  be  assumed,  and  if  the  cross 
sections,  (1,  are  known,  the  integrals  in  Equations  149  and  150 
can  be  evaluated,  at  least  numerically,  ahead  of  time. 


The  rate  of  loss  of  species  a  because  of  this  two- 
body  collision  is 


'  n(v  V  Vga^d“y  d“6  d“e  d“e 


(151) 


Consideration  of  the  integrals  in  Equations  149  through 
151  maket  it  clear  that  the  two-body  source  rates  must  satisfy  the 
following  "Compatibility  Relation" 


f  (2) 

JRae  d“a 

r  (2) 

■  *  • 

r  (2) 

Jr  _  du 

Y,ap  Y 

/<2)  _ 

r  (2) 

=  A,asdu6  *  . 

/R.  „Qdu 
/  e,ap  e 

Reduction  to  Spherically  Symmetric  Intermolecular  Potential 


That  the  above  definition  of  two-body  source  rates  is 

a  satisfactory  generalization  of  the  definition  given  in  standard  texts 

on  nonequilibrium  kinetic  theory  may  be  seen  as  follows.  Suppose 

that  no  nuclear  or  chemical  process  takes  place  during  the  two-body 

interaction;  then  particle  y  is  the  same  particle  as  particle  a, 

particle  6  is  the  same  particle  as  particle  0 ,  particle  e  does 

not  exist,  and  the  interaction  degenerates  into  a  two-body  scatteung. 

Because  the  integration  over  du^  reduces  to  a  delta  function  at  u  =0 

e  ’ 

Equations  149  and  150  give 
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(153) 


(2) 

RV  Kr7‘)  fffijpa?  *) ce(%- *■  *)  '  ga6 

‘^VV^pK^e^ 

If  the  intermolecular  potential  is  spherically  symmetric, 

n(Wgae)du6  =  n(“<  SaB)  d«)  (154) 

where  III  is  the  angle  of  scattering  (this  expression  is  merely  the 
statement  that  conservation  of  angular  momentum  in  an  elastic  two- 
body  collision  shows  that  the  scattering  cross  section  can  only 
depend  upon  the  relative  velocity  of  the  two  approaching  particles 
and  upon  the  angle  of  scattering  (see,  for  example.  Section  3.  3  of 
Reference  18).  Furthermore,  if  u,  and  ^  are  specified,  the  con- 
servation  of  momentum  allows  calculation  of  u.  (mathematically, 

H  =  0  except  for  the  curve  ua  =  U^.Cj).  Thus,  E,uations  15J 

and  154  become 


{Z)! 

\  M*)  • 

-  •few  0  ("•■««)  d«  as 


V  r*  7 


(155) 


In  a  similar  fashion. 


° (V  U6  *  Je :  gag)  d\  ^6  d"g  =  0  (w>  ga(3)  d u»  du  (1 56) 


Thus,  Equation  151  becomes 


~(2)  ff  h 

Rd  =  fd(Ua*  r*  7  *  l) 8ag  ’  ^  (do  8ag)  dw  du 

(157) 

For  an  elastic  two-body  collision,  specification  of  u 
causes  ua  and  ug  to  become  definite  functions  of  .  Thus,  EqJa- 
tion  155  can  be  written 
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(2) 

Ry  =  #aL5a(vV"’)'?'t]-SaS-nf’-Sasj 

■  SlVV  V“’)'7't]  d“6  ddl  <158> 

Now  interchange  the  indices  in  Equation  158  and  note  that  conserva¬ 
tion  of  momentum  requires  that  e  =  a 

6ag 

(2) 

Ra  (vM  = 

■  £6  [“6(v  >  “)• 7'  *]  dS  d1’  (1 59) 

Now  let  Y  “*  a?  to  represent  particle  a  after  the  collision,  and 
6  “*3' 5  Equations  157  and  159  give 


(2)  _ (2) 

R  -  R 
a  a 


ffla  V  ■  fa£e)  eae  n  (“•  sae)  ^  <16°> 


Equation  160  ,  except  for  minor  differences  of  notation,  is  identical 
with  the  source  rate  equations  usually  given  in  texts  on  nonequilibrium 
kinetic  theory  (see,  for  example.  Equation  3.  34  of  Reference  18). 
Thus,  Equations  150  and  151  are  satisfactory  generalizations  of  the 
source  rate  terms  in  the  Boltzmann  equation. 

c*  Source  Rates  for  Three-Body  Collisions 

The  three-body  source  rates  can  be  treated  in  a  similar 
fashion.  Thus,  assume  that  particles  a,  g,  and  y  approach  a  col¬ 
lision  region  with  velocities  ua.  up,  and  u  ,  and  define 


ma  +  m0u_  +  m  u 
q  q  g  g  y  y 
m  +  m.  +  m 

a  8  y 


(161) 


Let  the  cross  section  for  the  reaction 


a*  8;  Y  -*  6,  e 


(162) 


(163) 


where  the  notation  denotes  that  particle  Y  may  not  enter  into  a  chemi¬ 
cal  or  nuclear  reaction,  but  serves  only  to  carry  off  excess  energy. 

It  should  also  be  noted  that  particle  e  need  not  be  present  (u  =  6(0)) ; 
for  example,  in  the  case  of  molecular  recombination,  atoms  a  and  0 
recombine  into  molecule  6  with  particle  y  serving  to  conserve 
energy^  The  flux  of  particles  of  species  a  towards  the  collision  area 

18  fa(  V  r'  7  I  ua  "  u3  I  '  with  similar  expressions  for  the  fluxes  of 

particles  of  species  0  and  y  .  Then  the  rate  of  formation  of  particles 

6  and  e  due  to  this  particular  three-body  interaction,  T  is 

6e,a0y 

r6e,a0y(v  Vr#t)  =  -/Z/l  V*3I  *  I  V“3I  '  IV“3^ 

'  £a(v  **  *)  *  fp(V  *)  *  fyPY'  *) 

'  n(v'UvV'5 y)d\d*$dS  (164) 

The  rate  of  formation  of  any  one  species,  6  ,  by  three-body  processes 
is  then 

(3)  , 

Ra.a0;y(V*'V  =  /F  SC.aPyfvV1'4)^  <l65> 

The  rate  of  loss  of  particles  of  species  a  due  to  this  three-body  pro¬ 
cess  is 


•  n(v  \  I V  V“Y)d3Y  d36  d\  dSe  (166> 

“  iS  °Uar  fr°m  E<JuatlonB  164  through  166  that  the  three-body  source 
rates  must  satisfy  the  following  "Compatibility  Relation:" 
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/1(3>  -  /1<3>  /•(  3) 

J  a@,Yduo  =  /VyS  =  7*6,016! y^t 


■  / 


R  a  du 
e,ag;y  e 


(167) 


With  the  source  rates  in  Equation  144  defined.  Equa¬ 
tion  143  can  be  used  to  generate  the  conservation  relations  needed 

for  the  numerical  study  of  high-altitude  coupling  effects. 
d*  Equation  of  Change 

Let  xc7u)  rePresent  some  parameter  associated  with 
species  a.  Then  the  equation  of  change  is 


f  df  F  df 

I  Cl  -4  a  Cl  rr 

/xa  T7”  +  U  '  — qr  + - 

J  CL  dt  9r  m^  3u 


/  oata<ffi  *  A(£a*a)d3  * 


d"  =  A(V*a)d“  (168» 


=  £  ("a  <  Xa  >)  069) 


Similarly 


f  ^  hfa 

A au"aTdu  *  /g'  (faV)l!  *  *•(»,(](,»,))  (170) 

O 

f  aTr'ii  d"  =  f*5  (x^  £a)dS  'A  ^  '  MAa)® 

(171) 

fa)d“  ■  f,  ^  £a  _  maX  =  0  (172, 

*-  Jucx  . 

min 
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because 


fa(%riax)  "  fa('^amin)  =  0 

Also,  as  indicated  in  Equation  138 


(173) 


■A  •  F 
du  a 


(174) 


thus.  Equation  171  becomes 


/’  F  3f  f  f 

Xa~  '7§  d5  *  *  k—  •  V-v  d5 

a  9u  y  a  ma  u 


a  ,  ^  -•  . 

—  <F  •  V-.  y  > 
ma  a  u 


(175) 


Finally 


(176) 


*  ifx,]  (176) 

vhere  V|.Xa]  denotes  the  time  rate  of  change  of  X(l  per  unit  volume 
a"  (r* t}  because  of  chemical,  molecular,  and  nuclear  interactions. 

Using  Equations  168  ,  1 70  ,  175  ,  and  176  , 

s(VV)”'(WV)  *  ^<fa‘W  +  ’[xa]  <»» 

3-  CONSERVATION  EQUATIONS 

The  "Conservation"  equations  required  for  the  numeri¬ 
cal  work  are  obtained  by  choosing  ^  as  the  mass,  momentum,  and 

energy  of  a  particle  of  species  a.  To  handle  the  excited  states 

of  the  atomic  and  molecular  species,  it  is  convenient  to  enumerate  the 
possible  internal  energy  states  of  species  a;  thus,  if  state  a 

(B  =  1.  2,  .  .  .  ,  Na)  denotes  one  of  the  possible  internal  energy  states 
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of  species  p,  the  corresponding  internal  energy  is  denoted  by  ea 
(internal  energy  per  particle).  Furthermore,  state  itself  will  ^ 
have  a  distribution  function,  f^u,?  t),  which  gives  the  number  of 
particles  of  species  a  in  internal  energy  state  C(Xp  per  unit  velocity 
volume  and  unit  physical  volume  at  time  t.  It  is  clear  that  the  over¬ 
all  distribution  function  for  species  a  is  given  by 

Na 

*a(u,  r,t)  =  fap(u,r,t)  (178) 

Let  i|r  (u)  be  an  arbitrary  function  of  velocity  only. 


Then,  as  in  Equation  1  25 


<*a> 

/*(3>fada  ,  f 

'  /fadS  ‘  Tj*  fa du 

(179) 

Similarly 

<*v 

P 

(180) 

where 

% 

=  ff^Sl 

(181) 

Substituting  Equation  178  into  Equation  179 

f  Na 

"a<*a>  *  p 


Nn 

if 

B=ly 


f<xQ  du 


N, 


”aR  4cu> 

g=l  e  p 


(182) 


by  use  of  Equation  180 


By  use  of 


ma  =  m*3  - 

Equation  182  gives 


where 

%  =  m,s 

a*  Conservation  of  Mas. 

Let 

*a  =  ma  = 

Then 

■  0 

and  Equation  177  gives 

M(%ma)+7,  (”a(,IV'V) 

Summing  Equation  188  over  ^11  (X 


Multiply  Equation  182  by 

N 
a 

O  ( iir  N  -  V  ^  v 


: 


Equations  189  and  190  show  that 


dp 

lf+5'(p  a<V)  *  4[ma]  (191) 

where  A  [mj  is  the  time  rate  of  change  of  mass  of  species  a  per 

unit  volume  at  (r,t)  because  of  chemical,  molecular,  and  nuclear 
processes. 


b-  Conservation  of  Momentum 
Let 


Xa, 


(3 


rn  u' 
a 


J 


(192) 


where  u^  is  the  x?  component  of  u. 


Now, 


<*a  • 


”a<Fa-V> 

du 


ma<Fa6i> 

ma<Fa> 


Using  Equation  193  in  Equation  177 

£K<v) + 7  ■  (pap<%s>) 


(193) 


PaB  - 

—  <F0>  +  A[ma^e]  (194) 


Summing  Equation  194  overall  g  and  using  Equation  190 

+  7’(pa<V‘a>)  =  ^  < Fa}  +  7  [ma"a]  <195> 

where  A^uJ  is  the  time  rate  of  change  of  momentum  of  species 
a  per  unit  volume  by  chemical,  molecular,  and  nuclear  proces 


ses. 
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(200) 


+  -?■  >1 +  7 


■  [ftx^e^au  +  -r|)j 


p , 


Op 


=  ^<Fa‘S>+T4%+-r 


Summing  over  all  g  and  using  Equation  1 90  to  effect 
2  V* 

Pa/Ua  ^  ~  ^  *  Equation  200  gives 


4afa+(r))J+-(-T<Vcf>)^-§ 


PaB*V'V 


where 


P  e 

a  a 


But 


2 


where 


P 

m 


,a<f°'Sa>+B?lm“(eilB  +  -r) 


(201) 


N/ 


'ac 


fs{  Pas 


(202) 


^PVV1^  >  =  Va(V  + 


H. 


LPapeag<^p)-Paea<ua; 


=  paea<ua>  +  Z-Qa 


3=1 


Qac 


(203) 


(204) 


-g  r  pa0eap<ua3  - «a> 

is  the  excitation  heat  flux  vector  of  species  a.  Qa^  arises  from  the 
fact  that  the  excited  states  of  species  a  may  have  different  velocity 
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distributions  than  the  ground  state;  indeed,  because  different  excited 
states  result  from  different  sets  of  collisions,  it  seems  probable  that 
the  distribution  functions  of  the  various  excited  states  will  differ  dur¬ 


ing  some  nonequilibrium  period. 
'  V  =  0  and  Qig  =  o. 


U  %  =  IT  fa  *  then 

a 

In  general,  however,  a  diffusive 


heat  transfer  within  species  a  due  to  the  existence  of  excited  states 
with  unequal  distribution  functions  will  exist. 


Using  Equations  203  and  204  in  Equation  201  the 
energy  equation  becomes 


£(Va)  +  ’  >a< 


N„ 


-  V 


where 


(206) 
2\ 


s°  E  «a+(r) 

is  the  total  energy  per  unit  mass  of  species  a,  and  *  jm^  + 

is  the  time  rate  of  change  of  the  total  energy  of  species  a  per  unit 
volume  at  (r,  t)  because  of  chemical,  molecular,  and  nuclear  processes. 


Equations  191,  195,  and  205  are  the  required 
conservation  equations.  They  can  be  written  in  more  conventional 
form  by  use  of  the  pressure  tensor  and  the  heat  flux  vector.  The 
pressure  tensor  for  species  a  is  defined  as 

paBJ  = 


(207) 


at (pa <V)  +  (Pa  +  pa(  V <u^),  j  =  (Fa‘ >  +  A [maua] 

Ct 


Now  consider  the  term 


from  Equation  205 


(212) 


(213) 


The  heat  flux  vector  for  species  a  is  defined 


Pn 

a  /  1  J 

TF”  (  V  * 


~<T«V 'iV'ES 


(214) 


1  1 

vrt  v a. 
a  a  “j 


Thus 


(ua  -<ua>)  K  -  <4})  (^j  '  ^ 

(ui‘<uo>)(uo  -2ui<"Oj>+<“a>2) 


<vavav“j>  "  <uaua  >  "  *<»a^><«oj>+  <ua‘>  <ua>2 

-  <u^>  <ua2>  +  2<uj>  <u^>  (ua.)  .  <uo‘>  <ua> 

=  <ul(a«  ■  <\>))  - 

Using  Equation  213  in  Equation  213  (21 


(215) 


<v^v  ■  -2^-2<v^-<^>)('‘i-<ui>)> 


Noting  that 


(Ua-<ua>)  * 


-2<u»j><Uaj>(ua-<ua>> 


<ua>-<ud>  *  0 


(216) 


(217) 
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and  using  Equation  208  ,  Equation  216  gives 

<vavajva>  : 


T1  P1J 

-2  -  2<«  )  — 
aJ 


Substituting  Equation  218  in  Equation  214 


Na 


Q_ 


-  h  - 


'  -  pa  J<^j> +  Z 

J  B=1  B 


or 


*  -Qa-paj<"aj>  +  gS 

Using  Equations  213  and  219  in  Equation  205  gives 

^aea)+kea<“ai>).i  *  ‘  Qa.i  *  (paS>).i 

+  ^<FaS>  +  4[ma(a  +  :r)J 


Now  consider  the  term 
le  = 


<FaSi> 

N L 


■  /s 


.  faQ  F 1  u.  du 

rj  e  a  i 


(218) 


(219) 


(220) 


(221) 


As  discussed  immediately  prior  to  Equation  138  ,  the  only  velocity 

dependent  force  is  perpendicular  to  the  velocity,  so  Equation  221 
becomes 

,  iN" 

F 

1 


.  /*1Na 

C  /S  ^ B  u 

a  ^  p=l  P  1 

^  U-  du 
CL  I  CL  1 


Fa  <  Uai> 


(222) 
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becomes 


Using  Equations  221  and  222  ,  the  energy  equation 


(223) 


=  "  9  '  Sa  '  <%>  +  4  [ma^a  +  (223) 

where  4j%^a  +  -f-J  is  the  time  rate  of  change  of  energy  o'  species 
a  per  unit  volume  by  chemical,  molecular,  and  nuclear  processes. 


where  Ajm 


The  quantity 


wd  fo>> 

* 


Six  ■  %  j^v5’ 


(224) 


gives  the  rate  at  which  particles  of  species  a  carrying  property  cp^ 
are  generated  per  unit  volume.  Let 


reference  rate  of  particle  generatirn  per  unit 


volume 


(225) 


Then 


wa(l) 


■  2/1 

6=1  J 


s»  • 


(226) 


and  Equation  191  becomes 


.  +  ^  "  p  ( u  ) 

St  a ' 


KmaV1) 


(227) 


Similarly,  Equation  195  becomes 

+  5  ■  k<V  <\> +  pj 


pa<f?a/ma)+  Kmawa(S) 


(228) 


APPENDIX  II 
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CV{I,  J,  K) 

CV(I,  J.4N-3) 

CV(I,  J,  4N-2) 

CV{I,  J,  4N-1) 

DALPHA  (I) 

DEE  ZEE  (I) 

DELT 


ICMX;  1  sJs  MAXJC;  l  £  K  s  4  •  NSPEC) 

total  a-momontum  of  species  N  in  cell  J  in 
column  I  «PNV 'Wi+*J+*) 

total  r-momentum  of  specie.s  N  in  cell  J  in 
coiumn  I  «PNVuNr,i+i  j+4, 

total  energy  of  species  N  in  cell  J  in 
column  1  «PNVeN).+i.+i, 

sine  of  the  angle  between  rays  I  and  1+1 
0  s  I  s  ICMX)j(sin(ai+j  -  a.)) 

distanpe  along  a-axis  from  intersection  of  ray  I 
to  intersection  of  ray  I  +  1  (l  s  i  <  ICMX)- 

(zi+i '  V 

time  step 


DELTAR 
DELTAZ 
DFRND  (I,  N) 

DNODE(L,  J) 
DNODE(l,  J) 

DNODE(2,  J) 

DNODE(3,  J) 

DNODE<4,  J) 


r-inprement  along  arc  currently  being  treated 

z-increment  along  arc  currently  being  treated 

distance  along  ray  I  from  reference  point  to 
the  node  formed  by  the  intersection  of  floating 
boundary  N  with  ray  I  (l  s  I  s  IBMX;  2  sNs4) 

(l£Es4;l«Js:  MAXj) 

distanqe  along  ray  I  from  reference  point  to  the 
node  formed  by  the  intersection  of  ray  I  and 
arc  J  in  column  I 

distance  along  ray  I  +  1  from  reference  point  to 
the  node  formed  by  the  intersection  of  ray  I  +  l 
and  arc  J  in  column  I 

distance  along  ray  1+1  from  reference  point  to 
the  node  formed  by  the  intersection  of  ray  I  +  1 
and  arc  J  in  column  I  +  1 

distance  along  ray  1+2  from  reference  point  to  the 
node  formed  by  the  intersection  of  ray  1+2  and  arc 
J  in  column  1+1 


DRBND(I) 

r-increment  along  last  arc  in  column  I 
(l  s  I  s  ICMX);  (Ar.  ,  .  ) 

1+s»  Jmax 

DSTBLY(L) 

(1  £  L  £  2) 

DSTBLY(l) 

minimum  distance  between  the  arc  boundaries 
of  a  cell 

DSTBLY(2) 

minimum  distance  between  the  ray  boundaries 
of  a  cell 

DST1 

maximum  possible  ray- segment  length 

DST2 

maximum  possible  arc- segment  length 

DVOL(J) 

volume  change  of  cell  J  in  current  column 
during  time -step  (l  sjs  MAXJ) 

DZBND(I) 

z-increment  along  last  arc  in  column  I 
(i  *  I  *  ICMX);  (Az  .  ) 

1  *’  Jmax 

EB(N,  L) 

(1  «  N  £  NSPEC;  1  s  L  s  2) 

EB(N,  1) 

internal  energy  per  unit  mass  of  species  N  on 
cell  boundary  currently  being  treated 

EB(N,  2) 

internal  energy  per  unit  mass  of  species  N  on 

first  arc  of  next  column  downstream 

EBO(N) 

internal  energy  of  species  N  per  unit  mass  on 
the  boundary  for  a  problem  with  constant  boundary 
conditions  (l  £  N  £  NSPEC) 

EINT(L,  J,  N) 

internal  energy  per  unit  mass  of  species  N  in 
cell  J  of  column  I+L-l(l£L£2; 

1  ^  J  s  MAXJC;  1  s  N  s  NSPEC) 

GAM(L,  J,  N) 

ratio  of  specific  heats  of  species  N  in  cell  J 
of  column  I  +  L  -  1  (l  s  L  £  2;  1  £  J  <;  MAXJC; 

1  £  N  s:  NSPEC) 

GAMMA  (N) 

constant  ratio  of  specific  heats  for  species  N 
treated  as  a  perfect  gas  (l  £  N  £  NSPEC) 

GASCST(N) 

constant  gas  constant  for  species  N  treated  as 
an  ideal  gas  (1  £  N  £  NSPEC) 

IBCD 

nimber  of  cycles  between  which  BCD  output 
is  desired 

IBMX 

total  number  of  rays 

ICMX 

total  number  of  columns 

IEXTM 

flag  denoting  time-dependent  boundary  conditions 

IECXTM 

Is  0  constant  boundary  conditions 

(S:  1  time -dependent  boundary  conditions 

IGEOM 

geometry  option  word 

IGEOM 

At  plane  geometry 
|2,  cylindrical  geometry 

3,  spherical  geometry 

IMESH 

flag  denoting  whether  or  not  the  mesh  moves 

IMESH 

0  fixed  mesh 
la  1  floating  mesh 

(Note:  A  fixed  multi-column  mesh  in  which  the 
cells  differ  in  geometry  from  column  to  column 
must  be  treated  as  a  floating  meah.  ) 

IS  LIP 

flag  denoting  whether  or  not  arcs  must  be  continuou 
at  rays 

ISLIP 

0  ,  arcs  are  continuous  across  rays 
la  1  ,  arcs  are  not  necessarily  continuous  across 
rays 

KTAPE 

tape  number  for  BCD  input  tape 

MAXJ 

maximum  number  of  arcs  permitted  in  any  column 
of  the  mesh 

MAXJC 

maximum  number  of  cells  permitted  in  any  column 
of  the  mesh 

MTAPE 

tape  number  for  binary  restart  input  tape 

MXFARC 

maximum  number  of  floating  boundaries  permitted 
in  any  column  of  the  mesh 
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NARC(N,  L) 

NDIMEN 

NDIMEN 

NFRARC(I,  N) 


NMOMEN 

NMOMEN 

NSPEC 

NTAPE 

NTOT{L) 

P(L,  J,  N) 

PB(N,  L) 
PB(N,  1) 

PB(N,  2) 

PBO(N) 

RGAS(L,  J,  N) 


arc  number  of  the  free  arc  in  column 

NAKC(l.l)  =  NARC(1,  2)  =  1  (2  s  N  4;  ! 


I  +  L  -1 
s  L  £  2) 


number  of  dimensions 

1»  one  column 

2,  more  than  one  column 


arc  number  of  the  arc  which  represents  floating 
boundary  N  in  column  I.  NFRARC(I,  N)  =  0 
if  floating  boundary  N  does  not  occur  in  column  I 
(1  ^  I  *  ICMX;  2SNS4) 

number  of  momentum  equations  being  treated 

!li  radial  momentum  only 
2,  radial  and  axial  momentum 

number  of  species  under  consideration 

tape  number  for  binary  restart  output  tape 

total  number  of  free  arcs  in  column  I  +  L  -  1 

(1  £  L  £  2) 

partial  pressure  of  species  N  in  cell  J  of 
column  I  +  L  -  1  (l  sLs2;  1  s  Js  MAXJC; 

1  £  N  £  NSPEC) 

(l  £  N  £  NSPEC;  1  £  L  £  2) 

partial  pressure  of  species  N  on  the  cell 
boundary  currently  being  treated 

partial  pressure  of  species  N  on  the  first  arc 
of  the  next  column  downstream 

partial  pressure  of  species  N  on  the  boundary 
for  a  problem  with  constant  boundary  conditions 
(1  £  N  £  NSPEC) 


gas  constant  of  species  N  in  cell  J  of  column 
HL-l(lsL£2;l£js  MAXJC;  1  £  N  £  NSPEC) 
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T 


RHO(L,  J,  N) 

RHOB(N,  L) 
RHOB(N,  1) 

RHOB(N,  2) 
RHOBO(N) 
SFLUX(L,  J,  K) 

L 


density  of  species  N  in  cell  J  of  column 

.  "  1  J1  s  L  s  2;  1  j  s  MAXJC; 

1  s  N  s  NSFEc);(pa) 

(^Ns  NSPEC;  1  s  L  s  2) 

density  of  species  N  on  the  cell  boundary 
currently  being  treated 

density  of  species  N  on  the  first  arc  of  the 
next  column  downstream 

density  of  species  N  on  the  boundary  for  a 
problem  with  constant  boundary  conditions 

fluxes  Of  mass,  momentum,  and  energy  across 
the  jurfacesofa  cell  (l  s  L  s  5;  1  s  J  s  MAXJC; 

J.  s  X  fi  4  •  NSPEC/ 

1,  fluxes  entering  cell  J  in  current  column 
across  arc  J 

2,  fluxes  entering  cell  J  in  current  column 
across  arc  J  +  1 

3,  flpxes  entering  cell  J  in  current  column 
across  upstream  ray 

4,  fluxes  entering  cell  J  in  current  column 
across  downstream  ray 

5,  fluxes  entering  cell  J  in  next  column 
across  downstream  ray  of  current  column 


4N-3,  flux  of  axial  momentum  of  species  N 
K  J  4N-2,  flux  of  radial  momentum  of  species  N 

4N-1,  flux  of  energy  of  species  N 
4N  ,  flux  of  mass  of  species  N 


Sine  of  the  angle  between  ray  I  and  the  negative 
z-axis  (Isis  IBMX);  (sin  a.) 

SORCE(K)  (l  s  K  s  4  •  NSPEC) 

SORCE(4N-3)  z  body-force  on  species  N  in  cell  being  treated 


SORCE(4N-2) 


SORCE(4N-l) 

SORCE(4N) 

SVOL(L) 

T 

TBCD 

TBIN 

TIBIN 

TULT 
UB(  J) 

UFREE(L,  N) 
UN(L,  J,  N) 

UNB(N,  L) 
UNB(N.l) 

UNB(N,  2) 
UNBO(N) 


r  body-force  plus  radial  pressure  balance  on 
species  N  in  cell  being  treated 

work  and  heat  added  to  species  N  in  cell  being 
treated 

mass  source  of  species  N  due  to  chemical 
reactions 

volume  swept  out  by  arc  J  +  2  -  L  bounding 
cell  J  in  current  column  during  the  time- step 

time  (t) 

time  counter  for  BCD  output 

time  counter  for  binary  restart  output 

time  increment  between  which  binary  restart 
output  is  desired 

time  desired  at  end  of  run 

velocity  normal  to  itself  of  arc  J  in  the  current 
column  (l  s  J  s  MAXj) 

velocity  of  the  free  boundary  on  arc  N  of  column 
I  +  L  -  1  normal  to  itself  (l  s  L  s  2;  1  s  N  s  MAXJ) 

if  applied  to  an  arc:  velocity  normal  to  arc  J  of 
species  N  in  cell  J  +  L  -  2  of  the  current  column; 
if  applied  to  a  ray:  velocity  normal  to  ray  I  of 
species  N  in  cell  J  of  column  I  +  L  -  2 
(l  s  L  s  2;  1  s  J  s  MAXJ;  1  £  N  s  NSPEC) 

(l  s  N  s  NSPEC;  1  *  L  s  2) 

component  normal  to  boundary  of  velocity  of  species 
N  on  the  boundary  currently  being  treated. 

component  normal  to  boundary  of  velocity  of  species 
N  on  first  arc  of  next  column  downstream 


component  normal  to  boundary  of  velocity  of 
species  N  on  the  boundary  for  a  problem  with 
constant  boundary  conditions  (l  £  N  £  NSPEC) 
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UR(L,  J,  N) 

USTBLY(L) 

USTBLY(l) 

USTBLY(2) 
UT(L,  J,  N) 

UTB(N,  L) 
UTB(N,  1) 

UTB(N,  2) 
UTBO(N) 

UZ(L,  J,  N) 

VNODE(L,  J) 
VNODE(l,  J) 

VNODE(2,  J) 
VNODE(3,  J) 

1 


radial  component  of  velocity  of  species  N  in 
cell  J  of  column  I+L-l  (l  s.  L  s  2- 
1  s  *  MAXJ;  1  £  N  £  NSPEC) 

(l  sLs2) 

maximum  speed  of  a  characteristic  normal  to  an 
arc 


maximum  speed  of  a  characteristic  normal  to 
a  ray 

if  applied  to  an  arc:  velocity  tangential  to  arc  J 
of  species  N  in  ceil  J  +  L  -  2  of  the  current 
column; 

if  applied  to  a  ray:  velocity  tangential  to  ray  I 
of  species  N  in  cell  J  of  column  I  +  L  -  2 
(1  s  L  i  2;  1  «;  J  s  MAXJ;  1  sNs  NSPEC) 

(l  $Ns  NSPEC;  1  s  L  s  2) 

component  tangential  to  boundary  of  velocit  j  of 
species  N  on  the  boundary  currently  being  treated 

component  tangential  to  boundary  of  velocity  of 
species  N  on  first  arc  of  next  column  downstream 

component  tangential  to  boundary  of  velocity  of 
species  N  on  the  boundary  for  a  problem  with 
constant  boundary  conditions  (l  £  N  £  NSPEC) 

axial  component  of  velocity  of  species  N  in 
cell  J  of  column  I+L-l  (l  s  l  s  2; 

1  s  J  s:  MAXJ;  1  £  N  £  NSPEC) 

(l  £  L  £  3;  1  £  J  s  MAXJ) 

velocity  along  ray  I  of  the  node  formed  by  the 
intersection  of  ray  I  and  arc  J  in  column  I 

velocity  along  ray  I  +  1  of  the  node  formed  by  the 
intersection  of  ray  I  +  1  and  arc  J  in  column  I 

velocity  along  ray  I  +  1  of  the  node  formed  by  the 
intersection  of  ray  I  +  1  and  arc  J  in  column  I  +  1 
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Set  time  counters: 

BCD  output  counter:  TBCO  =  U. 
Binary  restart  output  counter:  TBIN 


The  MAIN  program 

sets  the  counter  indices,  calls  for  input,  calls  for  tin 
advancement  af  the  cells,  controls  the 
output  af  dota,  computes  the  time,  and 
competes  the  stable  time-step  far  the  next  cycle. 


Input  all  data  necessary  ta  run 
problem:  CALL  INPUT 


Set  the  time-step  parameters: 
STBL(I)  =  DST1 
STBL(2)  =  DST2 
USTBLYd)  =  0.0 
USTBLY(2)  =  0.0 


»  Begin  sweeping  through  mesh: 

DO  1  1C  =  1 ,  ICMX 

- 1 - 

Compute  the  time  required  far  the 

fastest  characteristic  ta  cross 
between  the  closest  arcs: 

TAUI  =  STBL(I  /USTBLY0I 

ta  advance  the  cells  in  column  1C: 
CALL  COLUMN 


Begin  search  for  the  minimum 
distances  between  the  rays,  STB L(2) , 
and  between  the  arcs,  STBL(I),  in 
the  mesh:  DO  2  ID  =  I,  N DIMEN 


No 


h 


STBL  (ID)  >DSTBLY  (ID)?  > 

Alt 


OSTBLY  (ID) 
is  not  a 
new  minimum 


The  minimum  distance 
in  the  current  column 
is  lass  than  the  minimum 
distance  in  any 
previous  column 


Set  the  new  overall  minimum: 
JTBL  (ID)  =  DSTBLY  (ID) 

I 


loop  not 
completed 


»|  2  CONTINUE  J— 
Loop  £  completed 

-|  1  CONTINUE  ~| 
^T^^^^ompleted 


Loop  not 
completed 


Another  cycle  is  completed.  Advance  j 

time  counters:  T  =  T  +  DELT;  J 

TBCD  =  TBCD  +  DELT;  TBIN  =  TBIN  +  DELI 


^  IEXTM  >0? 
There  are 
time-dependent 
Soundary  conditions 


Yes 


No  _ 

Boundary  conditions 
do  not  change 
with  time 


Update  the  unsteady  boundory 
conditions:  CALL  EXTIME 


{^TBCD  *  TIBCD?  V 

Nol  - 


Yes 
BCD  output 
is  desired 


Artput  BCD  data: 
CALL  BCDOUT 
TBCD  =  0. 


TBIN  s  TIBIN? 


Na 


k  Yes 

Out  binary  data: 

r  Binary  output 
is  desired 

TBIN  =  0. 

T  TULT? 


Yes 


No 


Computation 
is  finished 


rcAu’N 

.BOD; 


j - X - 

(What  is  N DIMEN? 


2-D.  Must 
also  consider 
ray  characteristic 


l-D.  .a  ray 
character  istic 


Compute  the  time  required  far  the 
fastest  characteristic  ta  cross 
between  the  closest  rays: 

TAU2  *  STBL(2)/USTBLY(2) 

Compute  the  two-dimensional 
stable  time-step: 

TAUI  *  (TAUI  e  TAU2)/(TAUI  +  TAU2) 


No 


limit  on  the  YES 

time-step 

Unsteady  boundary 
conditions  may 
limit  the  time-step 

Calculate  any  limit  which  the 
boundary  conditions  place  an 
the  time-step: 

CALL  TIMSTP 

_ i _ 

Set  the  time-sfe 
Reset  the  time-s 
STBL(I)  =  DSTI 
STBL<2)  =  DST2 
USTBLY(I)  =  0. 
USTBLY(2)  =  0. 

p:  DELT  =  TAUI 
ep  parameters: 

0 

3 

Figure  10.  Flow  Chart  of  MAIN  Program 
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<  j  <  jc?  Vs _ 

.>■"  r  L«t  oil 

Anofh*f  Yts 

, —  cell  »x?an< 


For  eoch  ipeclei,  compute  the  preiture,  deniity,  and 
normol  ond  tongentiol  component!  of  velocity  on  ore  J- 
DO  II  NS  =  I,  NSPEC 

CALL  WALL  (> '  I ,  J ,  NS),  A(1 ,  J I ,  MS),  P(I,J,  NS),  P(1  J|  NS1  URMit 
UNO,  J, .  ns:  .UNO.  J I ,  NS),  UTO ,  J I ,  Nsj,UT(2,  J I ,  NSJGAMfl  J  NS) 


For  ooch  ipeclei,  compute 
the  fluxes  of  man,  momenta 
and  energy  entering  the  lost 
r.ell  of  the  column  ocrou  the 
lost  ore: 

CALL  BCJAMX(IC) 


Compute  the  internol  energy  per  unit  men  for  each 
ipeclei  o,  ore  4:  CALL  THERM2(J) 


Compute  the  fluxei  of  mau,  momentum,  ond  energy 
of  each  ipeclei  entering  cell  0C,J)  oerpt.  ore  Jit 
CALL  FLUX  (2, 4,41, 1C,  AREARC(1,JI),  U6(JI)  I) 


Compute  the  eoureei  of  moit,  momentum,  ond  energy 
of  each  rpeclei  within  cell  (1C, 4)  during  the 
timo-itop:  CALL  SOURCE  (J) 


Begin  odvoncing  the  eelli.  Loop  through  the  ipeclei, 
40  5  NS  =  I.  NSPEC 


Loop  through  the  comervotion  equotloni  for 
ipeclei  NS:  00  4  NCE  *  |,4 


Compute  the  number  of  the  equotion  being 
odvonced:  NEQ  -  4  i  (NS.  I)  +  NCE 


^  What  li  NMQMEN7 
Rodiol  momentum  I  1 

°"ry  f 

OssiLC 

Eq.  NEQ  represents  Yes 
on  oxiol  momentum 
equotton  ond  is 
not  needs'  J  for 
this  colculofion 


Compute  the  net  rett  of  peneroMon  of  property 

Rodiol  ond 
oxiol  momentum 

No 

DO  9  NS1  •  1,  N5 

S^RCE(NEQ)  -  SORCE(NEQ)  *  SFLUX(NSI,  J,  NEQ) 

9  CONTINUE 

Eq.  NEQ  exists . 
Advonce  lr. 

j - 1 

1 

Advonce  property  NEQ  in  cell  (1C,  J): 

CV(IC,J,NEQ)  »  CV(IC,4,NEQ)  +  DELT  .  SORCE(NEQ) 

•  oop  nor 
compUtao 

loop  not 
completed 


6  CONTINUE  1+. 

Loop  |  compl »,  *e* 

S  CONTINUE  I 


_ Vw  >  *  I 

Utf  coll 

in  column  No  • 

>i  another 

,  ,  coif  In  ilWt  column 
Rmet  the  flue  Indlcei  for  tho  next  ueiii 

SFLUX(),4I,NEQ)  «  -  SFLUX(2.4.NEQI 


I  CONTINUE 


Loop  not  r  1  - R^RETURNj 

completed  com  .lefed  s.  _^/ 

Figure  11.  Flow  Chart  of  Subrout  ne  ADVNCE 


I  Lorp 
com  iletad 
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A*C<tC) 
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Subroutine  ARCGEO  computes 
tho  length  of  o  given  ore,  the 
line  and  cosine  of  the  angle 
which  the  arc  makes  with  the 
positive  Z  axis,  and  the  area 
af  the  cell  face  farmed  fay  the 


ARCGEOIIC,  \ 
JA.l,D0l,0D2)J 


id  •  ic  + 1  ■ 

JAMXC  »  JAMX(IC) 


Wh*  |,  NOIMfNT 


Compute  the  radial  and  axial 
Increments  af  the  current  arc 
from  those  af  the  previous 
arc  and  from  the  incremental 
distances  between  the  arcs 
clang  the  bounding  rays  of 
the  column: 

DELTAR  ■=  DELTAR  +  DD2  •  SN(ICI) 
-  DDI  e  SN(IC) 

DELTAZ  -  DELTAZ  +  DDI  .  CS(IC) 
-DD2  .  CSflCI) 


Rodlel  and  axial  Increments  af  arc 
given  by  INRUT  data: 

DELTAR  .  DR  END  (1C) 

DELTAZ  -  DZINDQC) 


Compute  the  length  pf  the  arc: 

ARCL(L.JA)  -  SQRTFfDELTAR  .  DELTAR  +  DELTAZ  .  DELTAZ) 
Compute  the  sine  and  the  cosine  of  the  angle 

which  the  arc  makes  with  the  positive  i  axis: 

ARCSN(l,JA)  -  DEITAVARCL(L,JA) 

ARCC5(l,JA)  *  D{LTA2/ARCL(L,JA) 


Plane  geometry  factor 
Is  unity: 

SI  •  1,0 


Cylindrical  geometry 


Compute  cylindrical  factor: 

LI  •  2  •  L  -  I)  L2  ■  II  *  I 
51  -  3.MI5R27  •  (DNODE(LI,JA) 
•  SN(IC)  r  DNODE(L2,  JA)  • 


Compute  spherical  factor: 

SI  -  I2.JW37  •  DNODE(L,JA) 
•  DNODE(L,  JA) 


Compute  the  area  af  the  cell  face  which  this 
orq  form*: 

AREARC(L.JA)  -  SI  •  ARCl(l,JA) 


Figure  13.  Flow  Chart  of  Subroutine  ARCGEO 
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Compute  next  column  number: 
iCl  *  1C  +  1 

Determine  number  af  cells  in 
column  1C:  JC1  *  JAMXflC)  -  1 
Determine  number  af  cells  in 
column  iCl :  JC2  -  JAMXflCl)  -  1 

jE 


Subroutine  COLUMN  computet  the  geometry 
of  each  cell  in  the  column  and  the  velacitiei 
of  the  free  aret  ond  nodei  in  the  column, 
appliet  the  boundary  condition!  an  the 
ft  ret  ond  loot  rayt  af  the  meth,  and 
col  It  for  the  flux  computation!  ond  the 
advancement  af  each  cell . 


No 

Fixed 

meth 


<iMESH  >  0?>*J - / 

"0 .  thiNi 


t>p  looting 
f  meth 


What  it  N 01  MEN? 


it  the  1  it 
column 


Compute  the  velocity  af  each 
free  nude  an  the  lent  ray: 

N1  -  NTOT(1) 

DO  26  N  -  1,  N1 
N2  »  NARC(N,  1) 

SiNGi  >  SN(IBMX)  . 
ARCCS(1  ,N2)  +  CSflBMX) 

.  ARCSN(1,N2) 
VNODE(2,N2)  * 
UFREE(1,N2)/SING1 
26  CONTINUE 


2-0 .  Another 
i  column  exiiti 


IMESH  >  0? 


No 


Fix  id  | 
meth 


Compute  the  fluxet  af  each 
tpeclet  acratt  the  latt  ray, 
ond  determine  the  fceteit 
characteriitic  normal  ta  the 
latt  ray:  CALL  BCIBMX 


Floating  )  Yet 
methj 

Calculate  the  geometry 
af  all  af  the  cellt  in 
the  next  column: 

|  CALL  GMTRY(2,IC1,2)  l 

— c  1 

Compute  the  density,  rodiol  and 
axial  velocities,  and  internol 
energy  for  each  species  in 
each  cell  af  column  ICl .  Compute 
the  velocities  af  each  species 
in  eoch  cell  af  columns  1C  and 
iCl  normal  and  tangential  ta  s 

ray  ICl:  CALL  THERMO(2, JC2,IC1) | 


Compute  the  fluxes 
af  each  species 
entering  eoch 
cell  af  the  first 
column  across 
the  first  roy. 
Determine  the 
fastest  character' 
istic  normal  ta 
the  first  roy: 

CALL  B0 1  (JC1 ) 

{ 


Fixed! 
moth  I 


Floating 
meth 


I  Calculate  the 
geometry  af 
I  af  the  cellt 
in  the  firtt 
column:  CALL 
GMTRY(i,  i.NDiMEN) 


Compute  the  density,  radiol  and 
axial  velacitiei, and  Internal 
energy  far  each  tpeciet  in 
eoch  cell  af  the  firtt  column: 
CALL  THERMO(t,JCl,l) 


Na 


Compute  the  velacitiei 
af  all  nodei  and  arct 
in  the  column,  and  the 
fluxet  af  each  tpeclet 
acratt  each  arc  and 
ray  face  in  the  column; 
advance  each  cell 
af  the  column;  and 
reiet  all  parameter! 
needed  far  the 
calculation  af  the 
next  column: 

CALL  DIFFEQOC) 


IMESH  >  0? 


Fix. 

met! 


K* 


Yet 


Floating 

meth 


Begin  computation  af  the  free  ~ 
velacitiei  in  column  ICl  and  he 
free  node  velocltlei  an  ray  ICl: 

DO  17  N  -  I,  MXFARC 

T - - 

Determine  which  ores  in  columns  1C 
and  ICl  represent  free  boundary  N: 

N1  *  NFRAKC<IC,N);  N2  -  NFRARCflCl  ,N) 


Free  boundary 
N  do«s  not 
exist  in 
column  ICl 


d 


Compute  the  velocity  of  arc  ] 
N2  normal  ta  Ittelf: 

CALL  BCFJN(ICI,2,N2) 


Compute  geometric  factor  for 
are  N2:  SING2  •  SN(ICI)  . 
ARCCS(2,N2)  r  CSQC1)  »  ARCSN(2, 


Begin  computation  af  the  free 
arc  veloeitiei  In  column  1  and 
free  node  petition!  an  ray  1: 

DO  8  N  -  1,  MXFARC 

-  - 

Determine  which  arc  In  column  1 
represents  free  boundary  N: 

N2  «  NFRARC(KN) _ 


Free  boundar^ 
N  does  not 
exist  in  column  1 


Yes 


0?^ 

Arc  N2  does 
represent  a 
free  boundary 


Compute  the  velocity  af  arc  N2 
normal  to  Itself: 

CALL  6CFJN(1 ,1  ,N2) 

r  * 

Compute  the  velocity  af  node  N2 
along  roy  1: 

VNODE(l,N2)  *  UFREE(1  ,N2)/(SN(1) 

•  ARCCS(1 ,  N2)  4  CS(1)  »  ARCSNQ ,  N2))  I 


completed  I 


N1 


Compute  the  dlitance  along  roy  1 
from  the  reference  point  to  free 
node  N  at  the  end  af  the  tlme-ttep: 
DFR ND(1,N)  -  DFRND(1,N)  + 
VNODE fl  .N2)  »  PELT 

B  ONTINUE  i-igM.W1 
—  ■  ■  ■  I  completed 


N1  >  0? 


Free  boundary 
N  does  not 
exist  In 
either  column 


Loop 

.  £2LJ 

completed  | 


17  CONTINUE 

Loop  lCam  - 
*  pitted 


H 


No 


>  07  > 

EEL 


F. o 


Compute  geometric  foctor  fo' 
ore  Nl:  SINGI  >  SN{IC1)  , 
ARCCS(1,NI)  +  CS(IC1)  . 
ARCSN(1,N1) 


Free  boundary 
N  exiiti  only 
In  column  ICl 


N2  >  07 


Yet 


Compute  the  velocity  af 
free  mde  N  an  ray  ICl: 
VNODE(3,N2)  - 
UFREE(2,  N21/SI NG2 


Free  boundary 

Free  boundary  | 

N  exists  only 

N  exiiti  in 

In  column  iC  i 

both  columns 

1  Compute  the  velocity  af  free  node 

N  on  ray  ICl: 

VNODE(2, Nl)  . 

UFREE(l,NiysiNGl 

U 


E 


Compute  the  velocity  af 
free  node  N  an  ray  ICl: 

VNODE (2, Nl)  x  (ARCL(1 , N 1 )  . 
UFREE (2,  N2)/SING2  * 
ARCL(2,N2)  *  UFREc(l,Nl)/ 
SiNGll/(ARCL(l,Nl)  + 

ARCL(2,  N2)) 

VNODE(3,N2)  =  VNODE(2,Nl) 


Figure  14.  Flow  Chart  of  Subroutine  COLUMN 
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OlFFEQflC) 


®*8'n  computation  of  velocities  of  ol!  nodes 
In  the  column,  Stort  with  the  nodes 
botwoon  ore  I  and  th«  next  fro*  ore: 

Ml  ■  I  Determine  the  totol  number  of  free 
orei  In  rhe  column:  N2  *  NTOT(I) 

* 

Bejln  coleuloting  the  velocitlei  of  the  nodes 
between  successive  pairs  of  free  ores: 

DO  25  N  -  2,  N2 


- - -  »<>  ■  wvs  vvmpjrii  me 

velocities  of  oil  of  the  nodes  on 
both  royi  of  the  current  column, 
computes  the  fluxes  of  eoch  species 
crossing  oil  the  roy  ond  ore  foces 
of  the  column,  odvonces  eoch 
cell  In  the  column,  ond  resets 
the  parameters  needed  for  the 
coicuiotlon  of  the  next  column. 


Determine  the  ore  number  corresponding  to 
free  hour  dory  N:  N3  -  NARC(N,  1). 

Set  Innermost  ore  to  be  computed  this 
pom:  N4  -  N3  -  1 


Arc  N4  Is 
also  o  free  * 
boundary 


,N4  »  Nl?] 


Determine 
node  velocities 


compl,ftdl 


Reset  outermost  free 
ore  index:  Nl  •  N3 


25CONT1NU! 
oop  I  completed 


Compute  the  fluxes  of  ipau, 
momentum,  and  energy  of  eoch 
species  entering  the  cells  of 
column  1C  ocrou  ray  ICI . 
Determine  the  fastest 
characteristic  normal  to 
toy  ICI:  CALL  RAYflC) 


* — ~  <IC  •  ICMX? 

Tr«o» 

ray  ICI  Yes  Roy  ICI 

already 

treated 


Compute  the  velocity  of  eoch 
ore  in  the  column,  the  volume 
change  of  each  cell  during  the 
time-step,  the  fluxes  of  eoch 
species  entering  the  first  cell 
oeroM  the  first  ore,  ond  the 
fastest  characteristic  relotive 
toonyorc:  CALL  ARCflC) 


Determine  the  number  of  cells  between 
free  ores  Nl  ond  N3:  FAD  «  N3  -  Nl 


Determine  the  incremental  velocity  changes 
between  nodes  an  royi  1C  ond  ICI: 

DVI  -  (VNODE(l,Nl)  -  VNODE(l,N3))/FAD 
DV2  -  (VNODE(2,Nl)  -  VNODE(2,N3))/FAD 


Set  the  outermost  ore 
to  be  computed: 

N3  «  Nl  +  1 
r  l**""  '  1 
Compute  the  node  velocities  due 
to  ores  N5  through  NX, 

DO  28  N6  «  N.  N4 
N7  *  N4  -  I 

VNODE(l,N«).  VNODE(l,N7)-DVI  I 
Vt  JDE(2,N«)-VNCDE(2,N7)-DV2 
28  CONTINUE 


Compute  the  fluxes  of  eoch  species  crossing  oil 
of  the  arcs  In  the  column,  apply  the  boundory 
condition  at  the  last  ore  of  the  column,  ortd 
odvonce  eoch  species  In  eoch  cell  of  the 
column:  CALL  ADVNCFflC) 


CompuK  the  free  node  positions  on  roy  ICI  ot 
the  end  if  the  time-step,  rezone  the  mesh 
if  necessory,  ond  reset  the  porometers 
needed  for  the  calculation  of  the  next  column: 
CALL  RESET  (1C) 


Figure  15.  Flow  Chart  of  Subroutine  DIFFEQ 


101 


Compute  the  equation  number; 
NEQ  «  N2  •  NS 


Compute  the  flux: 

PSI(N2)  -  AI(NS)  •  I1(N2,NS)  ♦  A2(NS)  •  ?2(N2,NS) 
SFLUX(IFLX,.II.NEQ)  ■  SFLUX(IFLX..I1,NEQ)  -  PS1QM2) 


Dow  ret  re  am 
ray  boundary 


Store  flux  entering  cell  (1C  +  I,  J2)t 

SFLUXQFLX  ♦  l,J2,NEQj  ■  SFLUX(IFLX  +  l,22,NEQ)  +  FSi(N2) 


Subroutine  GMTRY  compute,  the  node  position  of  eo«h  nod,  formed 
by  theocc*  in  column  iC;  the  length*,  ore  oa  of  revolution,  and  angle* 
from  tKi  Z  oxlt  of  eoch  ore  In  the  column;  the  volume  ond  crounectionol 
oreo  of  each  cell  in  the  column;  ond  the  minimum  perpendiculor  cell 
height  ond  width  in  the  column, 


Begin  computation  of  volume,,  area,, 
length,,  and  angle,  for  eoch  cell  In 
the  column:  DO  4  J  »  I,  JAMXC 


ChooM  Innermoit  untreated  ore  ond 
next  ore  out:  JARCI  »  JAMXC  +  I  -  J; 
JARC2  »  JARCI  -  I 


Compute  the  area,,  length,,  and  ongl,,s 
anoclated  with  ore  JARCI: 

CALL  ARCGEO(IC,  JARCI, L,DONODE(I,  NT), 
DDNODE(2,  NT)) 


1-0,  ore 
length  not 
important 


^  What  ItlDlM?  - 

>■  ■  i  r 
2  2-0,  f  te*t ^forjni^>imum  are  length 
^  DSTBLY(2)  >  ARCL(L,  JARCI)?  XY*‘ 
DSTBLY(2)  No 
olreody  I, 
the  minimum 


Chonge  the 
•poclng  Index 
to  telecl  the 
node  (pacing 

beyond  JARCI: 
NT  «  NT  .  I 


Set  the  ttoblllty  test  length,  equoi  to  the 
moximum  roy  Mgment  between  ore,  DSTI,  and 
to  the  moximum  ore  length,  DST2,  a,  given 
in  INPUT:  DSTBLY(I)  -  DSTI,-  DSTILY(2)  »  DST2 


Inltlollie  parameter,:  NFC  •  I;  NT  «  1- 
NFP  -  I,  ICI  -  1C  +  1 


Begin  computation  of  node  spacing,  on  roy* 
1C  ond  ICI  due  to  ore  in  column  1C: 

DO  I  N  -  2,  MXFARC 


Identify  ore  which  coincide,  with  free 
boundory  N:  NFA  -  NFRARCflC,  N) 


No 


JARCI  I,  tne 
first  cell  in 
the  column 


There  I,  a  cell 
corresponding  to  JARC2 


Compute  the  volume  and  crou 
Mctlonal  oreo  of  cell  (1C,  JARC2): 

CALL  VOLUME  (2,IC, JARCI, L,ODNODE(I, NT), 
DDNODE(2,  NT),  CELVOLQ.,  JARC2).  AREAfl. ,  J  ARC2)) 
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Loop  not  completed 


Call*  (IC,Jl)ond  (ICI,J2) 
have  their  entire  ray  I  Cl 
boundaries  in  common. 

S«f  fractional  paramoton: 
JC2MIN  -  Jl;  JC2MAX  -  Jl; 
PART(J1,J1)  -  1.0 


login  computation  of  tho  fluxoi 
on  tho  roy  IC1  foco  of  coll  (IC1,J1): 
DO  15  J2  -  JC2MIN,  JC2MAX 


login  swoop  through  special: 
DO  16  NS  -  1,  NSPEC 


Compute  tho  flaw  properties  of 
species  NS  on  tho  ray  IC1  foco- 
CALL  WALL(AC1,J1,NS),  A(2,J2,NS), 
P(I,J1,NS),  P(2,J2,NS),  0.,  UN(1,J1,NS), 
UN(2,J2,NS),  UT(1,J1,NS),  UT(2,J2,NS), 
GAM(1,J1,NS),  GAM(2,J2,NS),  RHO(1,JI,NS), 
RHO(2,J2,NS),  NS,  1) 


Compute  tho  downstream  and  upstream 
characteristic  speeds  of  species 
NS  relative  to  the  roy  IC1  boundary: 

CALL  STILT Y(A(1,J1, NS),  UN(1,J1,NS), 
0.,  USTILY(2)k 
CALL  STILTY(A(2,J2,NS), 

-  UN(2, J2, NS),  0.,  USTILY(2» 


■CSXEr 


Cell  (IC1,J2)  Is  the  lost 
cell  in  column  IC1  needed 
to  match  the  roy  IC1  foce 
of  cell  (IC.JI): 

JC2MAX  =  J2. 

Store  the  amount  of  the 
ray  1C  I  face  af  cell 
(IC!,J2)  paired  so  for: 
TSIDEO  -  TSIDEN 


All  of  the  available  roy 
ICI  face  of  ceil  (IC,JI) 

Is  needed  to  motch  the 
roy  ICI  foce  of  cell 
0C1.J2): 

PART(J1,J2)  -  1.  -  SPART 


16 CONTINUE 

completed 


Loop  not 
completed 


Compute  the  thermodynamic  properties  of  each 
species  on  the  boundory:  CALL  THERM2CI2) 


Compute  the  fluxes  of  each  species  on  the 
boundory:  CALL  FLUX(4,J1  ,J2,IC, 
PART(J1,J2)  ♦  AREARY0.J1).  0.0.  1) 


15  CONTINUE  I 


Loop  I  completed 


Select  first  uncompleted  cell  in 
column  ICI:  JC2M1N  -  JC2MAX 


login  computation  of  the  fraction  of  the 
cell  foce  on  roy  ICI  of  cell  (1C, J 1)  which 
also  bounds  cell  (IC1,J2).  The  mnemonic 
for  this  fraction  Is  PART(J1,J2). 

DO  9J2  «  1,  JC2 
PART(J1,J2)  -  0.0 
9  CONTINUE 


Set  fraction  of  foce  of  cell  (IC,J1)  which 
if  already  paired:  SPART  «  0.0 
Select  first  uncompleted  cell  In  column  ICI- 
J2  »  JC2MIN 


Compute  area  which  con  be  covered 
using  all  of  previously  covered  oreo 
of  cell  (ICI,  J2)  plus  oil  nan-paired 
roysurfoce  from  cell  (IC,J1): 
TSIDEN  -  TSIDEO  ♦  (1 .  -  SPART) 
*  AREARY(1,J1) 


Yes  Not  all  of  avoi fable 
ray  face  oreo  ,if 
,  .cell  (IC,JI)  is  needed 


Compute  the  fraction  of  the  roy  foce  oreo  of 
cell  OC,  Jl)  needed  to  motch  the  roy  foce  oreo 
of  cell  (IC1,J2):  PART(J1,J2)  - 
(AREARY3.J2)  -  TSIDEOl/AREARY(1,JI) 


Set  already  paired  oreo  to  zero  for  next  cell 
In  column  ICI:  TSIDEO  *  0.0. 

Compute  new  olreody  paired  fraction  of  ray 

face  of  cell  QC.Jik  SPART  -  SPART  4  PART(J1,J2) 


JC2  >  J2?  ; 
Select  Yet 
next  cell 

J2  ■  J2  4  I 


Put  ony  extra 
roy  oreo  in 
coll  JC2 


Loop  not 
completed 


6  CONTINUE 


Loop 

completed 


Figure  18.  Flow  Chart  of  Subroutine  RAY 
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RESET  PC) 


Compute  the  number  of  the  dawnstreom 
roy:  ICI  =  1C  +  1 


Bogin  computation  of  tho  poll  flora  of  tho 
fro#  node'  an  roy  ICI  ot  tho  ond  of  tho 
time-step:  DO  I  N  »  I,  MXFARC 

=  r== " 

Dotormino  which  ore  In  column  1C 
coincides  with  froo  boundary  N: 

NFA  «  NFRARCflC.N) 


Subroutino  RESET  computes  tho 
froo  boundory  node  positions  on  tho 
dowratroom  roy  of  tho  column  of  tho 
ond  of  tho  time-stop,  colls  for 
rozoning  of  the  mesh  if  necessory, 
ond  resets  tho  parameters  needed  for 
the  calculation  o.  tho  next  column. 


Arc  NFA  does 
represent  o 
free  boundory 


Free  boundary  N 
does  not  exist 
In  column  1C 


1C  -  ICMX? 


Tost  for  ond 
carry  out  needed 
rozoning: 

CALL  REZONE 


For  ooch  species,  reset 
gas  parameters  In 
cell  flCl.JA): 

DO  15  NS  -  I.NSPEC 
A(I,JA,NS)  «  A(2,JA,NS) 
EINT(I,JA,NS)  »  EINT(2,JA,NS) 
P(),JA,NS)  -  P(2,JA, NS) 
RHO(l,JA,NS)  -  RHO(2,JA,NS) 
RGAS(I,JA,NS)  »  RGAS(2,JA,NS) 
UR(I,JA,NS).  UR(2,JA, NS) 
UZ(I,JA,NS)  -  UZ(2,JA, NS) 
GAM(I,JA,NS)  «  GAM (2, JA, NS) 
15  CONTINUE 


Reset  the  geometric 
parameters  in  cell  (ICI.JA): 
AREA(I.JA)  .  AREA(2, JA) 
CELVOL(I.JA)  .  CELVOL(2,JA) 
AREARY(I.JA)  -  AREARY(3,JA) 


Free  node  N  on  ray  ICI  .esults 
from  on  ore  In  column  1C.  Set 
index:  NP  •  2 


Compute  the  dlstonee  olong  roy 
ICI  from  the  reference  paint 
to  free  node  N  ot  the  end 
of  the  time -step: 

DFRNDflCl.N).  DFRNDflCl.N) 
VNODE  (NP,  NFA)  .  DELT 


1  CONTINUE 


completed 


Free  boundary 
N  mo y  exist 
In  column  ICI 


- - <  ISLIP  >  0? 

A, a  may  | 

became  ^  Arc>. 
discontinuous  i  ,  cont' 

*  V  1C  <  ICMX? 


Determine  number  of  ores  In 
next  column;  JAMXC  *  JAMXflCI) 


[  i  -leiM  lm  m  w.hm 


Reset  ore  parameters: 

ARCL(I.JA)  «  ARCL(2, JA) 
ARCSN(I.JA)-  ARCSN(2, JA) 
ARCCS(I.JA)  •  ARCCS(2, JA) 
AREARC(I, J  ;  •  AREARC(2,JA) 
DNODE(I.JA)  -  DNODE(3,JA) 
DNODE(2,JA)  -  DNODE(4,JA) 

1 10  CONTINUE  Ll00^'”* 


Loop  I  completed 


Reset  the  number 
of  free  ores: 

NTOT(I)  •  I.fOT(2). 

For  eoch  species, 
reset  the  column 
parameters: 

DO  li  NS  *  I,  NSPEC 
EB(NS,  I)  *  EB(NS,2) 
PB(NS,  I)  •  PB(NS, 2) 
RHOBfNS,  I)  •  RHOB(NS,2) 
UTB(NS,  I)  *  UTB(NS,2) 
UNB(N3, 1)  .  UNB(NS,2) 

16  CONTINUE 


ore  parameter*: 

Nl  -  NTOT(I) 
DOI4N.  I,  Nl 
NARC(N,I)  • 
NARC(N,2) 

N2  =  NARC(N,I) 
UFREE(I,N2)  - 
UFREE(2,N2) 
VNODE(l,N2)  - 
VNODE(2,N2) 

14  CONTINUE 


Figure  19.  Flow  Chart  of  Subroutine  RESET 
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SOURCE  CIA) 


Plana 

9*o  ma  try 


There  is  no 
centrifugal 
pressure 
correction 
in  plone 
geometr  : 

X  =  O.b 


Who*  IsIGEOM?  y. 

2 

Cylindrical 

geometry 


Set  the  cylindrical 
geometric  factor: 

X  =  6.2831854 


_3 _ 

Sphericol 

geometry 


Set  the 
sphericol 
geometric 
factor; 

X  -  8.0 


Compute  the  radial  momentum 
source  far  eoeh  species 
due  to  the  centrifugal 
p  re  Mere  correction: 

DO  5  NS  ■  1,  NSPEC 
NEQ  =  4  *  NS  -  2 
SORCE(NEQ)  =  X  »  AREA(1,JA) 
*  P(I,JA, NS) 

5  CONTINUE 


^RETURN^ 


Subroutine  SOURCE  computes  the  sources  of  mass, 
momentum,  ond  energy  for  eoch  species  within  the 
given  cell .  This  version  of  the  subroutine  is  for  o 
chemi colly  inert  mixture  of  gases  which  do  not  collide 
with  eoeh  other. 


Figure  21.  Flow  Chart  of  Subroutine  SOURCE 


f  ocAic  A 
IC.  JA,  l.  DOI,  1 
002.  VOl,  AACAI) j 


D«fin«  nee ded  orithmelic  function!' 

|  SCRVIF(A,()  .  A  .  B 
SCRV3F(A,B,  C)  -  SCRV1F(A,B)  .  B  .  (A  •  C) 
SC8V4F(A,g,C,D)  «  (A  «  (fl  *  C))  t  (C  *  0) 


Subroutine  VOLUME  compute  the 
volume  between  two  prescribed 
orct  or  between  two  successive 
positions  of  the  some  ore.  Under 
given  conditions,  the  subroutine 
olso  computes  the  cross-tectionol 
oreo  of  o  prescribed  cell  of 
revolution. 


Compute  plone  volume  term*- 
X  =  0.5 

Y-  SCRVIF(DD1,SN(IC))  t 
SCRV1F(DD2,SN(1C)) 


(Compute  column  indices:  1 C I  s  IC  +  I- 

LI  =  2  *  L  -  lj  L2  =-  LI  t  | 

Compute  common  geometric  terms1 

VI  -  0D2  •  ONODE(L1 ,  JA)  -t  001  t  (002  t  DNODEfL2  till 
V2  -  0NODE(L I , JA;  .  DNOOE(L2.JA>  °DE(L2' JA)) 

T 

2  |  Cylindricol 


I  Compute  cylindricol  volume  term*: 

X  *  1.0471976 

Y  *  SCR  V3F (DO I , SN(1C), 2 . 0  *  DNODEfLI  JA))  t 
VI  .  SN(IC)  .  SN(IC1)  *  SCRV3F0D2 , 
SN(|CI),2,0  *  DNODL(«.2,  JA)) 

2  -  SCR  V4F(V1,D  NODE  (L1,JA), DD1.V2)  •  SNOC)  + 
SCRV4f(V2,DNODE(L2,JA),DD2,V2)  «  SNflri) 


Compijt.  the  volume  of  cell  (IC,  JA  -  |). 

VOL  .  X  »  (V  .  DEEZEEQC)  *  Z  •  DALPMAOC)) 


Compute  sphericol  volume  terms- 
X  <  4.1887790 
Y  *  0.0 

Z  =  SCRV4F(VI, DNODEfLI,  JA), 
DDI.V2) 


Volume 
ond  oreo 


Vthot  ii  IGEOM?  >- 
2  fCyiindricol 


__3 _ 

Sphericol 


Compute  sphericol  oreo  terms: 

W  =  1.5707964 
X  =  0.0 

Y  -  (2.0  •  DNODEfLI, JA)  +  DDI)  .  DD1 

Z  >=  0.0 
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I 


Subroutine  WALL  computet 
the  pressure,  density,  ond 
normal  and  tongentiol 
components  of  velocity  of 
the  given  species  on  the 
prescribed  cell  boundary. 


Compute  the  velocity  of  the  characteristic 
leaving  the  boundary  In  the  upstream 
cell,  VI,  and  in  the  downstream  cell,  V2: 
VI  *  Itv-i  -  Al;  V2  *  <JN2  +  A2 


< 


VI  <  u? 


At  feait 

one  choracterl.tlc 
propogote. 
up*troom  from 
the  cell 
boundary 


C 


Vei 


No  . 


Ye. 


V?  >  U? 
Characteristic, 
propogote  in 
•och  direction 
from  the  cell 
boundary 


> 


No  characteristics 
propagate  upstreom 
from  the  cell 
boundary 


No 


All  informotion  cn  the 
boundary  comes  from 
upstream.  Use  upstreom 
properties: 

PI(NS,L)  -  PI 
RHOB(NS,L)  «  RHOl 
UN((NS,L)  «  UNI 
UTB(NS,  L)  -  UT1 


N-  characteristics 
popagate  downstream 
from  the  cell  boundary 


Information  on  the  boundary  comes  both  from 
upstreom  ond  downstream .  Use  the 
characteristic  solution  for  the  boundory 
properties:  ' 

TERMI  .  Al  •  RHOl;  TERM2  «  A2  •  RH02; 
TERM3  «  TERMI  4  TERM2 
PB(NS,L)  ■  (TERMI  »  P2  4  TERM 2  •  PI  4 
(TERMI  •  TERM2  *  (UNI  -  UN2))VTERM3 
UN«(fS,L)  -  ((PI  -  P2)  4  TERM 2  •  UN2 
4  if  Ml  »  UNI1/TERM3 


Compute  the  houndory 
density  using  the  upstream 
cell  properties: 

RHOB(NS,L)  •  RHONEG(Pl, 
RHOI,GAMI,GAM2, 
PB(NS,U) 


Yes 


All  information  on  the 
boundory  comet  from 
downstream.  Use 
downstreom  properties: 
PB(NS,L)  -  P2 
RHOB(NS,L)  -  RH02 
UNB(NS,L)  -  UN2 
UTB(NS, L)  -  UT2 


Moss  crossing 
the  boundory 
comet  from 
upstream 


^  UNB(NS,  L)  >  u?  ^ 


No 


Mom  crowing 
the  boundory 
com#,  from 
dowmtrtom 


Sot  tho  tongont 
on  tho  bounder 
tho  tongentiol 
tho  upitroom  co 
UTB(NS,L)  *  l 

lol  velocity 
equal  ta 
reladty  in 

II: 

ITI 

Compute  the  boundory 
domity  uling  th* 
downstream  coll  properties: 
RHOB(NS,L)  .  RHONEG(P2, 
R  HQ2/  GAM2,  GAM  1 , 
PB(N$,L)) 


[ - 

Sot  tho  tongon 
on  tho  bounder 
tho  tongentiol 
downstream  co 
UTB(NS,L)  - 

lol  velocity 
y  equol  to 
velocity  In  the 

1: 

JT2 

Figure  24.  Flow  Chart  of  Subroutine  WALL 
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APPENDIX  III 


EQUATIONS  OF  RADIATION  TRANSPORT 


i.  GENERAL  EQUATIONS 

The  equations  of  radiation  transport  may  be  derived 
in  the  same  manner  as  the  equations  of  hydrodynamics  (Appendix  I). 
Thus,  the  photon  distribution  function,  fRv#(  r*  ,  c*  ,  t#) ,  denotes 
the  number  density  of  photons  at  point  (r*  ,t+)  having  velocities  in  the 

increment  dc*  about  c*  and  frequencies  in  the  increment  dv^ 

about  v*  (the  subscript  asterisks  denote  quantities  with  physical 

dimensions).  The  number  density  of  photons  at  point  ( r  ,  t  )  is 
then  * '  * 


nRv*  =  /fRV*G*»c#,t#)dc#  (2 

where  the  integration  is  carried  out  over  all  velocity  space.  The 
total  number  density  of  photons  of  all  frequencies  is 


nR* 


(231) 


Since  the  energy  carried  by  a  photon  is  h^ ,  the  total 

radiative  energy  density  at  point  ( r  ,  t  )  is 

eR*  *  J7h*V*  "Rv* dv*  (232) 

In  equilibrium,  however,  it  can  be  shown  that  the  value  of  eR|>  is 
(Section  4  of  Chapter  XI,  Reference  2)  * 


eR*  =  4<J*T*/c* 


(233) 
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(240) 


/h*V*^*fRv*d  c* 


Rv 


4  a 


T  4/v 

#  */ 


R„ 


Definitions  230  ,  237  ,  and  240  involve  moments  of 
the  photon  distribution  function.  This  distribution  function  is  given 
by  a  Boltzmann-like  equation  similar  to  the  Boltzmann  equation  in 
Appendix  I;  since  photons  cannot  accelerate,  however,  the  form  of 
the  equation  for  photons  is  slightly  simpler.  As  derived  in  Reference 

2  (Section  3  of  Chapter  Xi),  the  equation  for  the  photon  distribution 
function  is 


dt.  fRv*  +  C*  ’  V*fRv* 


/^RVjA 

where  the  term  ^  3t*  /interactions  accounts  for  all  changes  in  the 
distribution  function  due  to  scattering,  absorption,  and  emission. 

Integrating  Equation  241  over  all  values  of  c 

#  ’ 

■£^£rv*  d°*  fRv,S*  d3»  -A^)  (242) 


/3£Rv*\ 

V  St,  )i 


interactions 


(241) 


where 


V  •  c 


interactions 

(243) 


because  r*  and  c*  are  independent  variables.  Interchanging  the 
differential  and  integral  operators  and  normalizing  according  to 
Equation  237  ,  Equation  242  gives 


'Rv 


+  c  V  •  Q 
dt  RV 


■  m 


(244) 


interactions 
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I 


Taking  the  first  moment  of  Equation  241  by  multiplying 


ky  c#  and  integrating  over 


c#  gives 


at/VRv^,  +  As  *RV, ■  (^Avvt^_ 


Multiplying  by  h^VR  ,  dividing  by  40,7*  and 
240  gives  * 


*  /  interactions 

(245) 

using  Equation 


J-  n  ,  3  /=*SV*fRVn<'S 

at  Rv  *  d — I — 7 - ■ 

*  4a*T*/vR* 


(dt*  qrv): 


I  (246) 

interactions 


The  radiative  pressure  tensor  is  defined  as  (Equation  2.  1  3  ,  Chapter 
XI,  Reference  2) 


A  =■  Av 


c*  c*  - 
Rv^  c  7"  dc* 


(247) 


where 


(248) 


Since  pressure  is  an  energy  densitv  • 

_  gy  aensity,  it  is  convenient  to  normalize 

PRv  in  the  same  fashion  as  eRv;  thus,  as  in  Equation  237  , 


Rv  ~4~j 

4ff*T*/c*% 

Introducing  Equation  249  into  Equation  246 


(249) 


•  PM 


interactions 


~TZ —  +  cV  •  p 
®t  RV 


■  cm 


interactions 


(250) 
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2.  SPECIFIC  TERMS 


a.  Interaction  Termj 


The  interaction  terms  on  the  right-hand  side  of  Equa¬ 
tions  244  and  250  represent  sources  of  radiative  energy  and  heat 
flux  due  to  absorption,  emission,  and  scattering.  These  terms  depend 
upon  the  properties  of  the  material  through  which  the  photons  are 
propagating,  and  can  be  expressed  in  terms  of  absorption,  emission, 
and  scattering  coefficients.  The  following  discussion  of  these  effects 
is  based  upon  Reference  19  (Sections  3  through  5). 

b.  Absorption 

The  volume  rate  at  which  species  a  absorbs  photons 

of  frequency  v  is  given  by  the  product  of  the  photon  flux,  c  fn 

'  *  RV^  » 

the  absorption  coefficient  per  unit  mass,  *  va  ;  and  the  density, 

Pcv 


(d  **  f Rv*)ab  s  Pa*H  va*  c*  fRv 


(251) 


Multiplying  by  ,  integrating  over  velocity  directions,  and  using 

Equations  230  and  237  gives 


\  ^  /abB 


u;  E  Pa,  HVa,  c,  eRv 


(252) 


Introducing  the  Bouguer  number. 


Ju  =  P*H*L* 


(253) 


where  K  >Jc  is  a  reference  absorption  coefficient.  Equation  252 
becomes 
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(254) 


-  c  B 


u 


eRv  £  Pa 
a 


K 

va 


Treating  the  heat  flux  vector  similarly. 


/h*v*c*f Rv^dc* 

U.  ~l  ~ 

*  4a*T*/x 


£p 


a 


Ok, 


'va*  * 


cBuQ 


R 


v  T  pa»va 
a 


(255) 


c*  Spontaneous  Emission 


The  rate  at  which  unit  mass  of  species  a  emits  energy 
at  frequency  V  in  the  direction  ct  is  given  by  the  emission  coefficient, 
j  VCU  •  Thus, 


_a 

at. 


(fR''*h*v*)i 


^  Pa*  J*  va* 


(256) 


rsp  em 

Integrating  Equation  256  overall  c  dividing  by  4a.T4/c  VR 
and  using  Equations  230  and  237 


fe), 


gives 

*  a  ^0L>!(jva^ 


sp  em 


U,  „  4, 

*  4a  T  /c.vD 
*  *  Ra 


(257) 


where 


■'va, 


/ 


J  dc 
Jva.  * 


in  the  net  rate  of  emission  of  radiant 
mass  of  species  a.  Defining 


(258) 


energy  at  frequency  v  per  unit 
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(259) 


'VCL, 


'va 


l**a*T*  / ' 


C*VR, 


and  using  Equation  253  ,  there  results 


m 


sp  em 


c  Bu  £  Pa  j 


va 


(260) 


^  ^  va#  *s  no*  a  ^uncti°n  of  the  direction  of  energy 
emission,  as  seems  reasonable  for  spontaneous  emission,  there  will 
be  no  net  direction  of  the  emitted  energy;  thus,  there  is  no  change  in 
the  heat-flux  vector  due  to  spontaneous  emission: 


(*) 


sp  em 


(261) 


d*  Induced  Emission 


The  rate  per  unit  volume  at  which  species  a  is  stimulated 
into  emitting  photons  of  frequency  v#  in  the  direction  c*  is  given  by 
the  product  of  the  photon  flux,  CjJ(fRv  ;  the  coefficient  of  induced 
emission  per  unit  mass,  J  ;  and  the  density,  p  ; 

T  CC.'  • 


p‘*f) 

\  dt*  /In  em 


?  Jva*  C*  fRv 

a  *  *  ■rN-  Jje 


(262) 


Equation  (262)  is  similar  to  Equation  251  for  absorption;  thus,  by 
analogy  with  Equations  254  and  255  , 


'  'In  em 


c  B  e_  Y* 

u  Rv  p  J  K 

a  'a  va  va 


(263) 


and 
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(264) 


(268) 


is  the  cross  section  for  scattering  a  photon  into  any  direction.  Inte 
grating  Equation  267  over  all  c*  and  using  Equation  240 


c*h*v*^^vai{svaiit/fRviitS  dS 
u*  4a*T*4/vR 


■cB«5Rv^aV 


va 


(269) 


Equation  269  shows  that  scattering  cannot  change  the 
direction  of  the  heat-flux  vector.  Thus,  the  radiative  energy  approach, 
ing  a  scattering  region  may  be  divided  into  pencils  of  radiation  coming 
from  many  different  directions;  Equation  269  shows  that  each  of  these 
pencils  will  be  attenuated  by  the  same  attenuation  constant, 

C  sva*va  *  Consequently,  the  ratio  of  the  number  of  photons 

in  any  one  pencil  to  the  number  in  any  other  pencil  will  be  unchanged 
as  a  result  of  the  scattering.  The  net  heat-flux  vector,  therefore, 
will  not  change  direction  as  a  result  of  the  scattering. 


f.  Final  Equations 


Combining  Equations  244  ,  254  ',  260  ,  263  ,  and 
265  gives  the  equation  for  the  radiative  energy  density: 


i  de_. 

1  Rv  -• 

- - - +  V  *  Q 

C  St  Rv 


Bu?'>af‘vaeRv(1-Jva)-jva]  (270> 


Similarly,  combining  Equations  250  ,  255  ,  261  , 
and  269  gives  the  equation  for  the  radiative  heat-flux  vector: 


B“QRv?',a’*va(l-Jva+svoJ  (271) 
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\ 

\ 


8*  Energy  Absorbed  by  One  Spec 


The  net  rate  at  which  radiant 


energy  is  absorbed  into 


“  e"ergy  °f  Specie8  a  is  <*«-»  “y  integrating  Equation 
*  260  ,  and  263  over  all  frequencies: 


abs  by  a 


'  Jva)  '  jva]dv  (272) 


as  written  in  Equation  272  ,  the  radiative  energy  density  is  normalised 
With  respect  to  radiative  parameters,  that  is,  the  normalisation  is  as 

shown  m  Equation  234  .  The  hydrodynamic  energy  density,  however. 

is  normalized  by  the  term  a  Tj 

P#U)j<  •  Changing  the  normalization  of 

Equation  272  accordingly. 


(*Jr) 

\ dt  / 

'  € 


abs  by  a 


where 


■a*T*4/c* 

#  U*  2  RV  K  RV  ^  "  Jva  ^  ~  ^vaj 

l°Bli/|^RvH,Rv(1  ”  Jva)~  jvaj  dv  ( 


40 

P*U*3 


(273) 


(274) 


is  the  Boltzmann  number. 
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APPENDIX  IV 
MAXWELL'S  EQUATIONS 


For  convenient  computational  application  to  a  wide 
variety  of  problems,  Maxwell's  equations  must  be  normalized  into  a 
form  which  can  be  related  to  a  similarly  wide  variety  of  initial  conditions. 

RATIONALIZED  MKS  UNITS 

\ 

Maxwell's  equations  in  rationalized  mks  units  are  given 
in  Table  1-1  of  Reference  20  as 


e  v  •  E 

o*  *  * 

— 

% 

(275) 

-♦  -♦ 

V  •  B 

*  * 

= 

0 

(276) 

V  X  E 

(277) 

*  * 

at* 

= 

\(]* +  c°*tC) 

(278) 

Equations  275  through  278  are  to  be  normalized  into  a 
form  convenient  for  studying  the  flow  of  ionized  gases.  For  this  purpose, 
the  reference  for  the  charge  density  is  taken  as  the  reference  particle 
density  multiplied  by  the  charge  of  the  electron,  and  the  reference 
electric  field  intensity  is  taken  as  the  field  intensity  which  this  reference 
charge  density  would  generate  if  divided  into  positive  and  negative  charges 
separated  by  a  Debye  length.  Since  the  Debye  length  represents  the 
distance  of  charge  separation  at  which  the  energy  density  of  the  electric 


122 


fltId  eqUaU  016  energy  density  of  the  thermal  motion,  i,  j.  felt  that 
this  ohoice  of  normalization  will  yield  a  nondimensional  electric  field 
of  order  unity  throughout  the  gas.  The  magnetic  flux-density  is  then 
normalized  by  equating  the  energy  density  of  the  reference  magnetic 
flux-density  to  that  of  the  reference  electric  field  intensity. 


The  reference  charge  density,  q  ,  is  thus 


(279) 


Assuming  that  the  length  scale  for  variations  in  the  electric  field  is 
the  Debye  length,  L ^  ,  Equation  275  gives 


(280) 


N*qe*LD* 


(281) 


The  Debye  length  can  be  estimated  by  equating  the 
energy  density  of  the  field  E*  to  the  thermal  energy  density  of  a 
gas  without  internal  structure  at  a  characteristic  temperature  T 


e  T*: 


C  E* 

o#  * 


T  N*k*T* 


(282) 


Substituting  Equation  281  into  Equation  282  and  solving  for  LD  gives 


3k  T  e 

*  >f  o# 

N  q  2 


(283) 
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Choosing 


No  U 

*  e*  * 


(289) 


as  the  reference  current  density,  Equation  278  become 

=  M  |  N  q  U  ] 


n  N  q  L 

-v  y  ...  *  e»  D» 


c  c 

Os's  * 


B  = 


N.  q  L  rr 

*  e#  D *u* 


3E  \ 

5  7 


or 


v  X  B  = 

IL 

+  J_  _9E_ 

(290) 

c 

c  3t 

2.  GAUSSIAN  UNITS 

Maxwell's 

equations  in  Gaussian  units  with  the  current 

density  measured 

in  esu 

are  given  in  Table  1-1 

of  Reference  20  as 

v  •  e 

*  ~ 

4v% 

(291) 

— ♦  — 4 

•  B  = 

V  Jjc 

0 

—4 

(292) 

v*  X  £ .  = 

V  JjC 

_1_ 

5B* 

c* 

9t* 

(293) 

0# 

■  +  -L  SF* 

c*  St  % 

(294) 

As  for  the  rationalized  mks  units,  the  reference  charge 
density  is  defined  as 


*R*  =  N*q 


e* 


(295) 


For  the  reference  electric  field  intensity,  Equati 


on  291  then  gives 
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4fN  q  L_ 

*  e*  D* 

e#  D* 


(296) 


if 


o* 


1 

4ff 


(2  )7) 


in  Gaussian  units. 

/ 

Equating  the  energy  density  of  the  reference  electric 
field  to  the  reference  thermal  energy  density  of  the  gas, 

2 


E 

8n 


T-  N  .k.T 
2  *  *  # 


Substituting  from  Equation'  296  and  solving  for  Ln 

* 

1/2 


(-3k*T») 

|ik.Vo,| 


(298) 


(299) 


Computing  the  reference  magnetic  flux-density  by  equating 
the  energy  lensity  of  the  reference  electric  field  to  that  of  the  reference 
magnetic  field, 


or 


8n 


B, 


E&_ 

8n 


E*(*Vo*  > 


1/2 


(300) 
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7 


if 


o* 


4l7 


(301) 


in  Gaussian  units. 

Applying  normalizations  295  and  296  to  Equation  291  , 
V  •  E  =  £q 


(302) 


Similarly,  Equation  292  becomes 
V  ‘  B  =  0 

Using  normalizations  296  and  300  in  Equation  293  gi 


(303) 


gives 


or 


-r-  X  EE. 
L*  * 


7  x  E 


_L  E  ?**L. 

C*  *  ^t 
-L  SB 

c  8t 


(304) 


Choosing 


N.q  u* 

#  e*  * 


(305) 


as  the  reference  current  density,  Equation  294  become! 


or 


T-  X  E,  B  = 

LJ  * 

* 


V  X  B 


_ +iiE 

c  cL*  St 


iL  +  _L  _se_ 

c  c  at 


(306) 
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138 


where 


(313) 


The  normalization  of  the  work  term,  Equati 
identical: 


•  °n  308  ,  is  obviously 


F  •  u 
a  a 


Vl  E'  \ 


(314) 


In  Gaussian  units,  the  force  is  (Table  1-1,  Referenc 


e  20) 


qa4^*+  c*  *  (315) 

Inserting  normalizations  2,6  #  ^  <  ^  ^  int0  Equatio„  J1#  givei 


N  q  L  E 
*4e* 

P*U* 


x5(fVo/2] 


(316) 


Using  Equation  313  plus  the  fact  that 


e 

o*  o* 


(317) 


in  Gaussian  units,  Equation  316  become* 


r  (J ’ M 


(318) 


Smce  Equation  308  is  valid  for  the  work  in  Gaussian  units 
as  rationalized  mks  units,  it  is  clear  that  Equation  314  is  the 
correct  normalized  form  of  the  work  term  for  Gaussian  units  as  well  as 
for  rationalized  mks  units. 
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APPENDIX  V 

INTEGRATION  OF  THE  EQUATIONS 


The  integral  equations  developed  in  Subsection  III-l 
have  two  general  forms;  these  are 


(319) 


XdVXd‘ft*  +?'  *  +  ’X"  +  iz] 


(320) 


Equations  319  and  320  are  to  be  integrated  over  a  cell  volume,  V  , 
and  over  a  time  step,  At  ,  in  plane,  cylindrical,  and  spherical 
coordinates. 


1.  EQUATION  319 


The  first  term  of  Equation  319  may  be  written  as 


where  3<p/d t  is  assumed  uniform  throughout  the  cell,  and  the  cell 
volume  is  a  function  of  time  because  the  mesh  moves. 
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If  the  cell  surface  moves  with  velocity  w,  the  rate  of 
cnange  of  cell  volume  over  a  small  increment  of  cell  surface  is 


A 

W-  v  dS 


(322) 


where  v  is  the  unit  outward  normal  to  the  cell  surface.  Assuming 

that  property  „  is  constant  on  each  element  of  the  cell  surface  during 

the  time  step,  the  rate  a,  which  property  „  changes  due  to  the  motion 
of  surface  element  dS  is 


4' if)  =  <pw-  i 

Integrating  Equation  323  , 


(323) 


<p  W-  v  dS 


Inserting  Equation  324  in  Equation  321  , 


f'nf 

Ar  Jkt 


dt^ 

3t 


-  I1  dt  f  <pw.  v 
J  */tn  /s 

t„  o 


(324) 


(325) 


Replacing  the  surface  integral  by  a  summation  over  the  faces  of  the 
cell,  and  assuming  that  the  properties  on  each  face  of  the  cell  remain 
constant  during  the  time  step,  Equation  325  becomes 


/dv/dt|?  • 

•V  -it 


At  E  <pWnS 
cell 
faces 


where  is  the  component  of  W  normal  to  the  cell  face. 


(326) 
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319  gives 


Using  Green's  theorem  on  the  second  term  of  Equation 


dt  V  •  0  = 


dt  /’f  . 

At  Js 


where 


-*  A 
0  •  P 


>/ds  *  i'Z  *nS  (327) 
cell 

faces 


(328) 


and  0  has  been  assumed  constant  on  each  cell  face  during  the  time 
step. 


The  last  term  in  Equation  319  is  integrated  under  the 
assumption  that  the  source  term,  ,  i8  constant  throughout  the 
cell  volume  for  the  entire  time  step: 


ldvIdtj1  *  *X*  vw 


(329) 


But  using  Equation  324  , 


V(0)  + 


dS  W  •  v 


=  V(0)  +  t  T  w  s 

r  n 

cell 

faces 

Substituting  Equation  330  into  Equation  329  , 


XdvX 


•  *1  [V<*o 


(330) 


cell 

faces 


(331) 
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Combining  Equations  326  ,  327  ,  and  331  ,  the 
integrated  form  of  Equation  319  is 


tev>. 


(<pV)  +  At 

Eo 


■  V.J 


(332) 


2.  EQUATION  320 


Equation  320  in  a  vector  equation  and,  because  of  the 
rota, ton  of  the  unit  vectors  in  curvilinear  coordinates.  Green's  theorem 
■s  no,  applicable  to  vector  integrands.  Consequently,  is  necessary  to 
eparate  Equation  320  into  its  scalar  components,  and  apply  Green's 
theorem  separately  to  each  component.  As  can  be  seen  from  the  sym¬ 
metry  of  the  computational  cells  illustrated  in  Figure  1,  i,  suffices  to 
consider  the  x  and  z  components  of  Equation  320  for  plane  and 
cylindrical  geometry,  and  it  suffices  to  consider  only  the  x  component 
tor  spherical  geometry.  Consideration  of  the  remaining  components 
will  not  yield  new  information  because  of  the  symmetries  involved. 

A 

■Let  l  denote  a  unit  vector  in  either  the  x  or  z 
direction  .  Then  the  required  components  of  Equation  320  are  given  by 

/dVXdtP'  t  +  5xi+«.32]  =  0  ,333, 

Noting  that  I,  is  a  constant  vector, 


iA.  M 
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(335) 


A  -  =t  -  A  3 

1  *  7  *  x  =  V  •  U  *  x) 

and 

*  *  ^  x  /,J  =  W  •  ?  X  £  +  7  *  w  X  i  =  7  .  (60  X  t)  (336) 

Applying  Green's  theorem  to  Equation,  336  gives 


V  X  CO 


A  -* 

P  •  W  X 


A 

•l 


(337) 


by  use  of  the  vector  identities.  Combining  Equations  333  through 
335  and  337  , 


XdvZd#<M,+^] 


J " 3S  J"  dt  V 


.A  -»  -*  A 

U  •  x  +  w  x  i)  =  0 


(338) 


Equation  338  cannot  be  approximated  in  the  same  man¬ 
ner  as  Equation  320  bee  .use  the  integrands  contain  geometric  factors 
which  are  not  constant  throughout  the  cells  in  curvilinear  coordinates. 

Thus,  the  z  and  x  components  of  Equation  338  must  be  considered 
separately. 

a.  z  Component  of  Equation  338 

As  discussed  above,  the  z  component  of  Equation  338 
is  required  only  for  the  plane  and  cylindrical  geometries.  As  shown 
in  Figure  1,  however,  the  angles  between  the  z  axis  and  the  normals 
to  the  cell  sides  are  constant  for  a  given  cell  regardless  of  the  value  of 
y  or  0.  Consequently,  the  integrands  in  the  z  component  of  Equation 
338  are  constant,  and  the  equation  may  be  approximated  analogously  to 
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Equation  319.  Using  Equation  325  the  c  u  i 

S  nation  the  symbol  r  to  represent 

either  x  in  plane  geometry  or  r  in  cylindrical  geometry,  and  the 

symbol  9  to  represent  either  y  in  plane  geometry  or  9  i„ 

cylindrical  geometry,  the  a  component  of  Equation  338  becomes 


«.2v, 


"■V  *■{&[(<. 


faces 


-  (X  V  +  v 

\  zz  z  *zr 


"r  +  Vr)]  s  - 


K  v 

2z  t 


(339) 


where 


9  (340) 

tor  the  ceiis  shown  in  Figure  1.  „ere.  nr  is  the  direction  cosine 

between  the  outward  normal  to  ft.  celi  face  and  the  r  (or  x)  axis, 

and  11  z  is  the  direction  cosine  between  the  outward  normal  to  the  cell 
face  and  the  z  axis. 

b*  r  Component  of  Equation  338 

In  the  cylindrical  and  spherical  geometries,  reference 

t°  Figure  1  shows  that  the  angle  between  the  x  axis  and  the  normal 

to  a  side  of  the  cell  is  not  constant  over  the  entire  side  of  the  cell. 

Consequently,  the  detailed  integrations  must  be  carried  out  for  these 
two  geometries. 

For  plane  geometry,  however,  these  geometric  factors 

are  constant  on  the  sides  of  the  cells  and  tho  int-a  *  ^ 

us’  and  the  integrated  equation  may 

be  wiitten  down  by  analogy  to  Equation  339  : 


1  35 


**A  =  WAH?n[M^h 

faces 

"(x**1'*  +  *XzVz  ■  V*)}  -  *2kV3  <341> 

In  carrying  out  the  integration  for  the  cylindrical  and 
spherical  geometries,  it  must  be  remembered  that  the  radius,  r,  is 
never  negative  while  the  x  coordinate  may  be  negative.  Consequently, 

these  integrations  must  be  carried  out  only  over  the  half  cell  for  which 
x  is  positive. 

For  the  cylindrical  geometry,  the  half  cell  of  integration 

is  given  by 

.1  s  f)  <  1 

2  y  s  2  (342) 

Thus,  Equation  338  becomes 

L dr  L  izhn rd9  IA*  (*' C0B  6  -  *e si" 9) 

+  ^cos9-^sineMl-^rIdt 

+  f  dA  f  dt+  f  dl  f  B(.(,)  d0  f  dtl 
A(0  =  -ff/ 2)  At  J  C  J-n/2  JLt  J 

[(Vr  +  *TeVe  +  X^z)'08  9  ”  (Xfjr^r  +  Ve 

+  X9z^)  8  +  (<Vg  -  Vz)COB  9  -  (Vz 

VrWz)s^n  8j  =  0  (  343) 
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where  /  d|  denotes  the  line  integral  around  the 
JC 

is  the  value  of  r  on  this  path,  and  A  is  the 
of  the  cell  as  shown  in  Figure  2. 


the  cell  at  fixed  0  , 


cross-sectional  area 


On  the  face  A(Q  =  n/ 2)  , 


On  the  face  A(Q  =  -  n / 2)  , 


On  the  faces  -  tt/2  <  Q  <  n/z  > 


~  constants 


(344) 


(345) 


(346) 


(347) 


(348) 


Using  Equations  344  through  348  in  Equation  343  and 
performing  the  Q  integration, 


/tdt(XdrX  da  2#+^r)-  f  . 


/  dA  xnn 

A(0=ir/2)  00 


‘  4=-q/2,dAX^  +  2/RU,dt(vr 


v  -  00  »  )  l  = 
:  z  e  zjf 


(349) 
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Noting  that 


l 


dV 


f  dr  f  dz  f  r  d@ 
Ar  Az  JO 


(350) 


and 


1“  •  tC‘ 


ft(i)  d0 


for  the  full  cell,  Equation  349  becomes 

XdV  t  &  +  ^)+  IdS  Idt  far  +  X 


V 
rz  z 


-  w  V  )  -  Zlt  \  dt  /dA  V 

9  it  Ja  *eQ 


=  0 


(351) 


(352) 


The  integrated  form  of  Equation  352  can  be  written 
by  analogy  to  Equation  339  : 


+  “  {R  [(*r  -  T  *ZrK 
faces 

-  (KSr  +  *  V*)]S  +  2ffAXgg  - 

(353) 

The  integration  for  spherical  geometry  is  handled 
similarly  to  that  for  cylindrical  geometry.  In  the  spherical  case,  the 
half  cell  of  integration  is  given  by 

-  7r/2  £  0  <:  ff/2  and  -  ff/2  <;  <p  s  ff/2  (354j 

In  this  coordinate  system,  the  x  component  of  Equation  338  becomes 
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w>c 


r  cos 


V  dd  ^  (fr  C°S  ^  cos  0  -  sin  0 


^(p  Sln  <P  cos  6^  +  i<2r  cos  0  cos  Q  -  J \  sin  Q 

-  1>  sin  (p  cos  0  +  /  dA  +  / 

v  J  L*4/q__/->\  J.  „ 


A(6=ir/2)  *A(0  =  .^72) 


A/ 2  A/ 2 

+  7  *1*/  ,  rlcos  <£>  d© 

-ff/2  -ff/2 

A/2  A/2  - 

+  /- */2r2  w/2r*  C°S  *  d9J-  [(Vr  +  XJB 

+  V*)”'  ”  008  9  '  (v,  +  We 

+  V«)sin  6 '  ( V'r +  We +  xw,v),in  ” cos  9 

+  (Ve  "  vj cos  °cos  6  -  (“r%  -  oy^)  Bin  e 
-  -  t*iri/gjsin  <p  nos  ejj  =  0  (355) 


Tj  is  the  radius  of  the  inner  face  of  the  cell  and  r  =  r 


is  the  radius  of  the  outer  face. 


On  the  face  A(Q  =  u/Z)  , 


2  1 


(356) 


(357) 
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On  the  face  A(0  =  -  it/ 2)  , 


(358) 


(359) 


w/i  the  faces  r.  =  constant  and  r  =  constant, 

•md 


(360) 


constant 


(361) 


Substituting  Equations  356  through  360  into  Equation 
355  and  carrying  out  the  angular  integrations. 


H C 


2  2h  r  j  r 

rdri_T^  i.  / 


dA  X, 


£<8=ir/2)  99 


-/ 


7V(0=  -ff/2) 


dA  X00  + 


[rl2t,xl 


+  r_  V 


2  r2 


Noting  that 


4*  </V2r2dr 


(362) 


(363) 


(364) 


for  the  full  cell,  Equation  362  becomes 


(365) 


The  integrated  form  of  Equation  365  can  be  written  by 
analogy  to  Equation  339  : 


»rv> 


facei 


“  X  V  I S  +  8  Av  -  Jr  v  > 

rr  r|  *00  ^r  tDj 


4 


(366) 


By  letting  r  and  0  represent  x  and  y  in  plane 
geometry,  defining 


(367) 


in  spherical  geometry,  and  defining 
V 


(0  in  plane  geometry 
rr  in  cylindrical  geometry 
4  in  spherical  geometry 


it  is  possible  to  combine  Equations  341  ,  353 


(368) 


and  366  into  one  form: 


faces 


-  *  Xr^z  - 


S  +  2”A  *99  - 


OJ 

(369) 
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appendix  VI 
CHARACTERISTICS 


The  field  properties  on  the  cell  faces  are  evaluated  by 
the  use  of  characteristics.  The  theory  of  characteristics  for  the 
hydrodynamic  equations  of  a  perfect  gas  in  equilibrium  is  given  in 
Reference  4  .  In  this  appendix,  the  theory  of  characteristics  is 
extended  to  a  mixture  of  gases  out  of  equilibrium,  to  the  equations  of 
radiation  transport,  and  to  Maxwell's  equations. 


1.  HYDRODYNAMIC  EQUATIONS 


The  gas  in  each  sell  of  the  mesh  is  assumed  to  be  at 
a  constant  state;  this  state  is  characterised  by  a  pressure  tensor 
defined  by  (see  Equation  207  ,  Appendix  I) 


P 

a 


V 


v  v  > 

0!  a 


(370) 


end  a  kinetic  temperature  defined  by  (see  Equation  2, 13,  Chapter  IX, 
Reference  2), 


3 

2 


k*T  * 
*  a* 


) 


or 


1_ 

2 


/  2 
m  (v 

a  a 


(371) 


The  normal  pressure  of  species  a  is  defined  by 


ot 


ip  ;f 

3  a 


a  t  -  V 

-T~  (v  •  v  > 
3  a  a 


Equations  371  and  372  give 


a 


n  K  T 
a  a 


(372; 


(373) 


as  the  kinetic  relation  for  the  species  ™  in  n  • 

6  BPecies  ot  in  the  cell  in  question. 

Now  consider  a  plane  with  a  constant  flow  on  each  side 
as  depicted  in  Figure  25 

A 


Figure  25. 

,  A  A 

where  n  and  s 

tial  to  the  plane. 


Surface  Coordinate  System 

are  unit  vectors  in  the  directions  normal  and  cangen 
respectively. 


In  terms  of  these  coordinates,  Equation  19  is 


tt+  ^(W>)  =  Vl) 

where  Equation  25  has  been  used,  and  Equation  20  is 

(W) +  d  k<v>2  +  POT„)  = 

..  t 

where  Equation  26  has  been  used. 


(374) 


(375) 
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But 


Pann  “  T  Pa  :  1  Pa  (376) 

where  the  assumption  of  constant  flow  within  the  cell  implies  that  the 
pressure  tensor  is  isotropic. 

'^ie  first  Law  of  Thermodynamics  gives 

%  *  wv(rj  (5^7) 

This  version  of  the  first  law  (with  the  various  species  decoupled)  is 
valid  for  the  flow  across  a  boundary  because  the  effects  of  collisions 
are  volume  effects  and  disappear  from  the  integrals  for  sufficiently 
small  cell  size.  Thus,  the  several  species  may  be  considered  to 
flow  across  the  cell  boundaries  without  interfering  with  one  another. 

Now 


(379) 


(380) 
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or 


de 


a 


cva  dTa  +  (a(l/pa)JT  d(1V 

O' 


c  dT  + 
va  a 


1*  1^2.1  .  P 
o  JI  I  a 
\  alp 
Ha 


d(l/, 


or 


c  dT  -  T  dS 
va  a  a  a 


'  alp 


a 


a 


Thus 


fe). 


a 


T  l*p  \ 

— SL  Ll 

c  ST 
va  \  a  Ip 


a 


By  use  of  Equation  373  , 


/SP 


a 


ST 


alp 


n  K 
a 


a 


Thus,  Equation  387  becomes 


ST 

w)s 


a 


T 

— -  n  X 
c  a 
va 


-  P  /c 
a  va 


Now  consider 

/  3Pa  I 


'  IWs 

v  a 


8  lpaXTa  1 

3<J/PaM  / 


a 


_ 

X 

T 

dp 

'a 

r  ST 

a 

m 

a 

X 

a 

»u/p0) 

»  „ 

+  P 

s  ° 

rt 

3(l/o  ) 

‘  OL 

*  m 

.a)  (385) 


(386) 


(387) 


(388) 


(389) 
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-t  p2 

or  a  c 


KT  o 

ot  a  i  + _ 


m  c 
ot  vet 


(390) 


Using  Equation  373  in  Equation  390  , 


_ 2L_  ,  , _ 

(i/p  )  \  me 

a  '  a :  vet 


y  P 

WT)  <391> 

Ot 


where 


m  c 
a  vet 


(392) 


Ignoring  volume  effects,  which  do  not  contribute  to  the 
surface  flux  terms.  Equation  374  becomes 


ap«.  ;  ,  ae 


+  <  v>  it  +pai;  <u«„>  = 0 


(393) 


If  the  changes  between  regions  (1)  and  (2)  are  small 

(ensured  by  sufficiently  small  cell  size!  the 

7  niaii  ceil  size),  the  flow  processes  will  be 

nearly  isentropic.  Thus, 


(394) 


8pa/9(1/V 

p7^v)Ja  ,395) 
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and 


(396) 


ad/pa> 

ap 

a 


a 


Thus,  Equations  391  ,  395  ,  and  396  give 

2 


dp, 


a 


3P 


a'  S 


a 


where 


a 


"P 


a 


y  P  p 
a  cr  a 


(y  P  /P  ) 
\  o  a  a/ 


1/2 


is  the  local  speed  of  sound  for  species 


a 


_1_ 

2 

i 

a 


(397) 


(398) 


a 


at 


Similarly 


Using  Equation  397  ,  Equation  396  becomes 


at 


(399) 


a 


dp, 


a 

an 


Defining 


M 

a 


ap 

_ a_ 

an 


(400) 


a 


<u  > 

an 

a 

a 


(401) 


Equation  393  becomes 
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~  +  a<v  <ua„>  +  =  o 


M  BP 
cl  a 


0a  an 


(402) 


In  the  same  way,  Equation  375  becomes 


M  dP  /  2  \ 

a  a  ^  3  ,  M  +  l  dP 

o  a  ~TT  +  TT  )  +  2<u  >  —  <u  )  +  I  _ £L  _  n 

P(X  a  9t  an  an  Bn  an  p  Bn  ~  ° 

'a 


Adding  and  subtracting  Equation  402  to  Equation 


(403) 


403  : 


a  M  ±  1  BP 

”  <a“n>  +  ~P^~  +  a«K  *  *)  £  <Uan> 


/m  +  1  ±  M  )  BP 

A-  a  a/  a 


a 


Bn 


=  0 


(404) 


Multiplying  Equation  402  by  ^  and  aubtracting  iron,  Equation 


404 


at  <Uan>  ±pa 3t  +  a«(Ma  ±1)^n{\n) 

(l  ±M  )  BP 
.  \ _ al  a 


P  Bn 

rCL  ° 


=  0 


(405) 


Rearranging 


at 


<uan>  ±  TV"  “IT  +  a  (M  ±  i)  f— <u  > 
a  a  3t  J  a '  a  /  [  3n  an 


l  3P  1 
1  a 


pa  Bn 
^a  a 


=  0 


(406) 
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N 


Define 


P  dP 


on* 


<u  >  if 

J  n  a 


P  a 

'a  a 


(407) 


Then  the  final  characteristic  equation 


is 


bJ 


+  a  (m  j 

ot  a  \  a  I  bn  a 


(408) 


Over  an  interval  of  constant  a^  and  M ^  ,  Equation  408  has  the 
solutions 


a± 


a± 


n  -  a  (m  ±  l\  t 

ot\  a  I 


Thus»  Ja+  is  constant  along  the  phase  paths 
n  “  aa  i^a  +  ljt  =  constant 


or 


"  =  no  +  aa(Ma+1)* 

and  Ja-  is  constant  along  the  trajectories 
n  "  aQ  “ljt  =  constant 


or 


■v.K-1)* 


(409) 


(410) 


(411) 


RADIATION  TRANSPORT  EQUATIONS 


The  cnaracteristics  for  the  equations  of  radiation 
transport  are  found  in  the  same  manner  as  for  the  equations  of 


hydrodynamics.  The  mathematical  derivat 
■n  the  present  case,  however,  because  the  < 


:on  appears  much  simpler 
equations  are  linear. 


As  in  the  case  of  the  hydrodynamic  equations,  the  volume 
source  terms  do  not  enter  into  the  determination  of  surface  properties, 
ns,  assuming  constant  radiative  properties  in  a  ceil  and  using  die 

same  (n,  s)  coordinates  as  used  for  the  hydrodynamic  equations, 
Equations  25  and  26  become 


=  0 


(412) 


c  aeR„ 

+  3  —  =  0 


(413) 


Equation  35  has  been  used  in  reducing  Equation  26  to  Equation  413 
Multiplying  Equation  412  by  a  constant  a 

_  7  constant.  a  ,  and  adding  and  subtracting 

Equation  41 3  to  the  resulting  equation  gives 


aert  *  Q 
i  Rv  WRpn 


)  +  c  aTT  (‘ 


aQ_  ±  -4- 
t  Rim  3 


(414) 


Rear 


ranging 


St  (eR„  ±  QRun  /a)  ±5 T  ({ 


e_  i  3a  Q 
^  R^  Rim, 


The  arguments  of  the  differential  operators  in  Equation  415 
identical  if 


)  =  0  (415) 


7^ 


(416) 
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Defining 


'Rt, 


n/T  Q 


Rim 


(417) 


and  using  Equations  416  and  417  in  Equation  415  gives 


9J 


V± 


an 


0 


Equation  413  has  the  solutions 

V  =  ct/yj) 

Thus,  is  constant  on  the  trajectories 


(418) 


(419) 


n  =  nQ  +  Ct/y/T 

and  is  constant  on  the  trajectories 


n 


nQ  -  Ct/y/T 


(420) 


(421) 


3*  MAXWELL'S  EQUATIONS 

As  in  the  case  of  the  hydrodynamic  equations,  the 
characteristics  for  Maxwell's  equations  are  found  by  neglecting  source 
effects  in  the  equations;  thus,  the  appropriate  forms  of  Equations  29 
and  30  are 


6B 

at 


+ 


c  7  X  E 


and 


0 


(422) 


aE 

at 


-  c  7  x  B 


0 


52 


(423) 


The  presence  of  the  vector  product  in  these  equations  necessitates 

the  defining  of  a  rectangular  Cartesian  coordinate  system  on  the  surface 
of  the  cell  as  depicted  in  Figure  26: 


cell  surface 


Figure  26.  Three-Dimensional  Surface  Coordinate  ; 
As  illustrated,  the  coordinates  (n,  .j,  .  )  form  a  right-handed 


System 


system. 


Remembering  *at  the  assumption  of  uniform  field  pro¬ 
perties  on  each  side  of  the  cell  surface  implies  that  the  only  derivative, 
which  exist  are  a/at  and  6/5n  ,  Equations  422  and  423  reduce  to 


(424) 


(425) 


(426) 


(427) 
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Adding  and  subtracting  Equations i  424  and 


Similarly,  Equations  425  and  426  give 

(Bs2  *  E*i)  *  C  IT  {\  ±  Esj  *  0 

Defining 


'l=fc 


B  ±  E 


Equati.cn  428  shows  that  JJ+  is  constant  along  the  path 


no  -  ct 


while  is  constant  along  the  path 


n0+  ct 


Similarly,  defining 


Equation  429 

n 

while  J2  is 
n 


Bs  ±ES 
2 


1 


shows  that  is  constant  along  the  path 


n0+  ct 


constant  along  the  path 


n0  '  ct 


427  gives 

(428) 

(429) 

(430) 

(431) 

(432) 

(433) 

(434) 

(435) 
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APPENDIX  VII 
GEOMETRICAL  FORMULAE 


The  formulae  for  the  cell  geometry  are  developed  as 
outlined  in  Subsection  III-5c.  To  minimize  loss  of  signifi¬ 
cant  figures,  all  distances  along  a  ray  are  measured  from  arc 

j  =  jmax  (the  clo8e8t  arc  to  the  z  axis).  It  is  convenient  to  define 
the  distance  between  two  successive  arcs  (  Figure  2  )• 

Ad  d  j 

i,j+l/2  i,j  ‘  di, jfi  (436) 

The  radial  and  axial  increments  on  arc  j  in  column  i+1/2  are 

iri+l/2,j  '  4ri+l/2,j+l  +  Adi+l.j+l/2  ’“‘“i+l 

‘  idi,j+l/2  8in<1i  (437) 

and 

izi+l/2,j  =  42i+l/2.j+l  +  4di,j+l/2  cosai 

'  “i+l.j  +  I/2  COsai+l  (438) 

The  length  of  arc  j  in  column  i+1/2  is 

iu/2,j  =  j(dri+1/2,j)2  +(izi+l/2,j)2  (439) 

and  the  angle  that  this  arc  makes  with  respect  to  the  positive  z  direc¬ 
tion  is  given  by 


convenient  to  define  the  following  auxiliary  form- 


=  idi.j+l/2  sinai 
T2  =  d,  .-r,  M . 


‘.J+>  i+I»  j+1/2  *  di+l,j+l\j+1/2 


+  Ad.  .  .  Ad 

i,  J+1/2  U”m,j+J/2 


(Vi  *  (Il>i^i.j+,^di,j+1/2).inai 

^  +  Vdi, j+1/2 


=  d.  .  ,  d 

l*J+l  i+1,  j+1 


(446] 


o  ,  1  <2  .iI;„  In  ter'13  °f  the6'  auxiliary  formulae,  the  voiume  of  ceil 
(1  +  1/2,  j  +  1/2)  is  given  by 


i+1/2,  j+1/2  ~  8j(VjAz  +  V^sin/ia) 

vhere  in  plane  geometry 


(447) 


1/2 

P,).  +  Pj.., 


(448) 


in  cylindrical  geometry 


gl 

=  tt/3 

(453) 

V! 

(454) 

V2 

=  (»4),  »ina.  +  (34)1+j  *taai+1 

(455) 

Az 

zi+l  "  zi 

(456) 

sin  Aa 

=  .in«x1+,  .  a.) 

(457) 

and  in  spherical  geometry 

h 

=  4n/3 

(458) 

vi 

=  0  ; 

(459) 

v2 

(460) 

Az 

=  0 

(461) 

sin  Aa 

=  1 

(462) 

The  surface  area  of  arc  j  in  column  i+1/2  is 

Si+l/2,j  '  Sl  S2 

(463) 

where  in  plane  geometry 

s,.  ■ 

1 

(464) 

S2  * 

(465) 

in  cylindrical  geometry 

si  ■ 

n(1i+l,j  8inai+l  +  di,j  sinai) 

(466) 

(468) 


and  in  spherical  geometry 


S,  =  4tt  df 


1 


S2  =  1 


J 


(469) 


The  area  of  the  face  of  cell  (i+1/2,  j+1/2)  on  ray  i  is 

Si,j+l/2  "  glAdi,j  +  i/2  (470) 

where  in  plane  geometry 


gl  =  1 


in  cylindrical  geometry 

BI  =  "(2di,j+l  +  4di,j+l/2),toai 
and  in  spherical  geometry 

g!  =  ° 

Note  that  spherical  geometry  does  not  require  S 

j  4 1  /  2 

The  cross -sectional  area  of  cell  (i+l/2,  j+1/2) 

Ai+l/2,j  +  l/2  *  g2  S3  +  63  S4 

where  in  plane  geometry 

g2  =  g3  =  0 

(plane  geometry  does  not  require  A  \ 

i+1/2,  j+1/2 

in  cylindrical  geometry 

V/2 


(471) 


(472) 


(473) 


is 


(474) 


(475) 


g2  =  (Zi+l  *  Z- 


g3  =  sin<ai+l  "  ai)/  2 


(476) 

(477) 
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APPENDIX  VIII 
STABILITY  ANALYSIS 


The  equations  of  radiation  transport  in  characteristic 
form,  Equations  418  of  Appendix  VI,  and  Maxwell's  equations  in 
characteristic  form,  Equations  428  W»d  429  of  Appendix  VI,  are  all 
of  the  general  form 


dt  3r 


0 


where 


c  >  0 


(483) 


(484) 


Integrating  Equation  483  qver  a  time  step  and  over  the 
volume  of  a  cell  in  a  fixed  one  "dimensional  mesh  gives 


n+1/2  s  Atfjn+*  rnl 

^ *  J*iHl  U  C  AtK*  "  J  J  (485] 

where  cell  (n+1/2)  is  bounded  by  cell  surfaces  n  and  n+1,  lowered 
indices  denote  properties  evaluated  at  the  beginning  of  the  time  step, 
and  raised  indices  denote  properties  evaluated  at  the  end  of  the  time 
step. 


To  study  the  growth  of  any  one  Fourier  component 
of  J±  ,  it  is  convenient  to  assume  an  exponential  variation: 

J±n+l/2=  J±«P» 


k(n+l/2)  Ar  -  wmAt 


t 


(486) 
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where  k  is  the  wave  number  and  u  i.  the  frequency  of  the  Fourier 
component  under  study. 


According  to  Equation  483  ,  J+  propagates  in  the 
positive  r  direction  while  J_  propagates  in  the  negative  r  direction. 
Thus,  the  characteristic  boundary  propertiei 


ss  are 


n 


and 


.n 


_n-l/2 


rn+ 1/2 


(487) 


(488) 


Substituting  486  through  488  into  Equation  485  yield,  the  two  equations 

.  ♦  0,6t(l  -  e  -**■)]  (48„ 


+  n+1 , 


and 


-  n+1/ 2 


_n+l/2 


f  -  c  6t(,ik4r  .  ijj 


(490) 


T 

a 

b 


Using  the  notation 

=  c  At 

=  1  -  cos  k  Ar 

=  sin  k  Ar 


(491) 

(492) 

(493) 


Equations  489  and  490  can  be  written  in  the 


matrix  form 


-jn+l/? 

+ 

»  _ 

^+  n+l/2 

= 

G 

• 

jn+i/2 

; 

J  .  '  / 

- 

- 

_  -  n+l/2 

(494) 


where  the  amplifies 


matrix,  G  ,  is  given  by 


1  +  r(a  +  ib) 


I 


1  +  T(a  -  ib) 


(495) 


The  matrix  G  is  a,  noymal  matrix,  since,  if  G*  denotes 
the  Hermitian  conjugate  of  G,  it  is  ofeviOHS  that 

GG  ••  G*G  *  0  (496) 

Since  G  is  a  normal  matrix,  the  necessary  end  sufficient  condition 
for  the  stability  of  Equations  485  is  that  the  eigenvalues  of  the  matrix 
G  satisfy  the  relation  (Equation:  4f  25  ,  Reference  14  ) 

I M  s  1  +  Of  At)  /497x 


or,  for  the  present  application 
XX  s  l 

The  eigenvalues  of  the  matrix  G  are 


(498) 


1  +  T(a  ±  ib) 

Consequently,  for  either  eigenvalue, 


XX 


* 


1  _ 

[1  +  r(a  +  ib)][l  +  r( a  -  ib)] 


(499) 


(500) 
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« 


(501) 


and  requirement  498  gives 

1  *  (1  +  aT)2  +  bV 

Because  a  a  0  according  to  Equation  492.  inequality  501  i5  satisfied  for 
a“  T  2  0  '  ThU!’  Eliaati°n  «5  is  stable  for  any  positive  ,ime  step. 
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